#!/usr/bin/env python3 import os import os.path import math import numpy as np # to install this lib, run as root: aptitude install python-numpy """ Novaco displacements for pentagonal symmetry, started from novaco_v0p7.py """ def abs1(x): return np.linalg.norm(x) def abs2(x): return np.dot(x, x) def cabs2(c): return c.real**2+c.imag**2 def rotation(angle): cangle = np.cos(angle) sangle = np.sin(angle) rotArray = [[cangle, -sangle], [sangle, cangle]] return rotArray def omegaL(q): return cL * abs1(q) def omegaT(q): return cT * abs1(q) def dyn_mat(q): sqrt3 = np.sqrt(3) qxa = q[0] * acol qya = q[1] * acol cosxa = np.cos(qxa) cosxa2 = np.cos(qxa/2) cosyas3 = np.cos(sqrt3*qya/2) offdiag = sqrt3*np.sin(qxa/2)*np.sin(sqrt3*qya/2) return k/M * np.array([[(3 -cosxa2*cosyas3 - 2*cosxa), offdiag], \ [offdiag, 3 - 3*cosxa2*cosyas3]]) def omega(q): #return eigenfrequencies and eigenvectors m = dyn_mat(q) w = np.linalg.eigh(m) eigenfreq=np.sqrt(abs(w[0])) # potentially dangerous, no failure if negative eigenvals eigenvectors=np.transpose(w[1]) sign = np.dot(eigenvectors[0], epsilonT(q)) sign = sign/abs(sign) eigenvectors[0] = eigenvectors[0]*sign sign = np.dot(eigenvectors[1], epsilonL(q)) sign = sign/abs(sign) eigenvectors[1] = eigenvectors[1]*sign # here it may help to do something about the eigenvectors phase continuity... return eigenfreq, eigenvectors def inBZ(q, rotatedtau): # returns True of False if q is in BZ, and a weight to make it # fuzzy around 0 inbz = True # the squared length of the b_i vectors /2 norm2tau_half = abs2(rotatedtau[0])*0.5 norm2_adjust = [norm2tau_half,norm2tau_half,norm2tau_half, norm2tau_half-1e-15,norm2tau_half-1e-15,norm2tau_half-1e-15] for i in range(len(rotatedtau)): # the projection of q in the rtau direction, times the rtau length aaa = np.dot(q, rotatedtau[i]) if aaa > norm2_adjust[i]: # violation of the 1BZ condition inbz = False cutoff = 1.0 absq = abs1(q) if absq < qmin: cutoff = math.pow(absq/qmin, 2) return inbz, cutoff def inBZ_dangerous(q, rotatedtau): # returns True of False if q is in BZ, and a weight to make it # fuzzy around 0 inbz = False rtaus = rotatedtau/2 rtaulen = abs1(rtaus[0]) for rtau in rtaus: inbz = inbz or abs(np.dot(q, rtau)/rtaulen) <= rtaulen cutoff = 1.0 absq = abs1(q) if absq < qmin: cutoff = math.pow(absq/qmin, 2) return inbz, cutoff def epsilonL(q): return q/abs1(q) def epsilonT(q): return [-q[1], q[0]]/abs1(q) def novaco_displacement(args): import re global acol, apot global pi, twopi global qmax, qmin global cT, cL global k global M pi = math.pi twopi = 2*pi rcatino = 1170. nsymm = 5 # default or commandline-modified values: apot = args.apot acol = args.acol alpha = args.alpha print("# apot:", apot," acol:", acol," alpha:", alpha,"degrees") # NOTE: LAMMPS harmonic potential is K_{lammps}(d_i - d_{i-1})^2 # in the notes the harmonic potential is K_{notes}/2 (d_i - d_{i-1})^2 # # K_{lammps} = K_{notes}/2 => K_{notes} = 2*K_{lammps} k = 2e-7 # in fkg / µs² (K_{notes}) k = k * 1e6 # in fkg / ms² M = 31.06 # in fkg cT = np.sqrt(3*k/M) * acol / math.pow(2, 3/2) # µm/ms cL = np.sqrt(k/M) * 3 * acol / math.pow(2, 3/2) # µm/ms V = 1e-2 # in zJ Vtrue = V/float(nsymm**2) # note the factor nsymm**2 in the proper definition of V! qmax = twopi/acol # *2./math.sqrt(3) qmin = twopi/rcatino print("# qmin, qmax:", qmin, qmax) klength = twopi/apot print("# klength:", klength) firstkappa = [klength, 0.] kvectors = [] for i in range(nsymm): gg = np.dot(rotation(twopi*i/nsymm), firstkappa) kvectors.append(gg) gvectors = [] for ki in kvectors: for kj in kvectors: gg = ki-kj if abs1(gg) > 0: gvectors.append(gg) if args.debug: print("debug1", gvectors) b1 = np.array([twopi/acol, twopi/acol/math.sqrt(3)]) b2 = np.array([0., 2*twopi/acol/(math.sqrt(3))]) maxlength = 3.1*klength maxradius = maxlength*math.sqrt(2) nmax = 7 tauvectorsu = [] for n1 in range(-nmax, nmax): for n2 in range(-nmax, nmax): tau = float(n1)*b1+float(n2)*b2 if args.debug: print("debug0", tau) if abs1(tau) < maxradius: tauvectorsu.append(tau) # first shell of reciprocal tau vectors, to be used in inBZ firstshelltau = [] for ishell in range(6): tau = np.dot(rotation(twopi*ishell/6.), b1) firstshelltau.append(tau) firstshelltau = np.array(firstshelltau) print("# particle-type x y ux uy |u| angle") # for ialpha in range(int(round(args.minalpha/args.deltaalpha)),int(round(1.+args.maxalpha/args.deltaalpha))): # en_1ph = 0. # alpha = ialpha*args.deltaalpha rot = rotation(-alpha*twopi/360) rotatedg = np.dot(gvectors, np.transpose(rot)) for filen in args.filenames: print("#now reading",filen,"one equil. position at a time") count=0 if filen=="-": f=sys.stdin else: f = open(filen, 'r') for line in f: line=re.sub("^\s+","",line).rstrip() nums=re.split('\s+',line) R0=np.array([float(nums[0]),float(nums[1])]) # the equilibrium position displac=np.zeros(2) # reset the displacement vector # print ("quiii", R0, "quaaa",displac) for gg in rotatedg: for tau in tauvectorsu: q = gg-tau inbz, cutoff = inBZ(q, firstshelltau) if inbz: omegas, epsilons = omega(q) for pol in range(2): scal=np.dot(gg, epsilons[pol])/omegas[pol]**2*np.sin(np.dot(q, R0)) displac+=scal*epsilons[pol] displac*=(-Vtrue/M) moddisp=abs1(displac) angdisp=np.arctan2(displac[1],displac[0]) print(1,R0[0],R0[1],displac[0],displac[1],moddisp,angdisp) if __name__ == "__main__": import sys import argparse commandname=sys.argv[0] desc="""Novaco energetics physics for pentagonal symmetry in k space OUTPUT: the angular dependence of the energy per particle by Nicola Manini & Giuseppe Santoro """ epil=""" v. 0.1 by Nicola Manini, 19/02/2024""" parser = argparse.ArgumentParser( formatter_class=argparse.ArgumentDefaultsHelpFormatter , description=desc, epilog=epil) parser.add_argument( 'filenames', nargs='*', default=['-'], help='Files with equilibrium coordinates. If none given, stdin is used') parser.add_argument( '-a', dest='acol', type=float, default=5.8, help='lattice spacing of particles (microm)') parser.add_argument( '-b', dest='apot', type=float, default=5.2, help='lattice spacing of corrugation (microm)') parser.add_argument( '-d', dest='debug', type=int, default=0, help='raise the debug level') parser.add_argument( '-g', dest='alpha', type=float, default=0., help='the mutual rotation angle (degrees)') parser.add_argument( '-i', dest='deltaalpha', type=float, default=0.1, help='UNUSED the angular increment (degrees)') parser.add_argument( '-z', dest='minalpha', type=float, default=0., help='UNUSED the initial angle (degrees)') parser.add_argument( '-m', dest='maxalpha', type=float, default=6., help='UNUSED the maximum angle (degrees)') parser.add_argument( '-x', dest='X0', type=float, default=0., help='relative x displacement (microm)') parser.add_argument( '-y', dest='Y0', type=float, default=0., help='relative y displacement (microm)') ## End arg parser definition args=parser.parse_args(sys.argv[1:]) d = vars(args) # adding prog to args (for unknown reasons it's not there already...) d['prog']=parser.prog # here the actual function doing the job is called: novaco_displacement(args)