path: root/nov_displacement_v0p1.py

#!/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)