forked from lijiext/lammps
337 lines
9.9 KiB
Python
337 lines
9.9 KiB
Python
#!/usr/local/bin/python-2.5/bin/python
|
|
|
|
Info="""
|
|
Module name: cfg2lammps.py
|
|
|
|
Author: (c) Andres Jaramillo-Botero
|
|
California Institute of Technology
|
|
ajaramil@wag.caltech.edu
|
|
Project: pEFF
|
|
Version: August 2009
|
|
|
|
Reads in an eff .cfg file and produces the corresponding lammps data and input files
|
|
|
|
NOTE: Unsupported functions will be reported in the output log
|
|
"""
|
|
|
|
# import essentials:
|
|
import sys, os
|
|
from math import log10
|
|
from shutil import rmtree
|
|
from getopt import gnu_getopt as getopt
|
|
import numpy
|
|
|
|
def printHelp():
|
|
print Info
|
|
print "Usage: python cfg2lammps cfgfile\n"
|
|
return
|
|
|
|
general="""
|
|
# Created %s
|
|
|
|
# General parameters
|
|
|
|
variable sname index %s
|
|
log ${sname}.log
|
|
|
|
units electron
|
|
newton on
|
|
boundary %s
|
|
|
|
atom_style hybrid charge electron
|
|
|
|
read_data data.${sname}
|
|
|
|
pair_style eff/cut %s
|
|
pair_coeff * *
|
|
|
|
thermo %s
|
|
thermo_style multi
|
|
"""
|
|
#%(date,name,boundary,cutoff,period)
|
|
|
|
minimize="""
|
|
# Minimization
|
|
|
|
min_style cg
|
|
dump 1 %s xyz %s ${sname}.min.xyz
|
|
dump 2 %s custom %s ${sname}.min.lammpstrj id type x y z radius fx fy fz rf
|
|
min_modify line quadratic
|
|
minimize 0 1.0e-5 %s %s
|
|
|
|
undump 1
|
|
undump 2
|
|
"""
|
|
#%(group,period,group,period,iterations,fcalls)
|
|
|
|
single_pt="""
|
|
# Single point energy
|
|
|
|
run 0
|
|
"""
|
|
dynamics="""
|
|
|
|
# %s Dynamics
|
|
|
|
timestep %s
|
|
|
|
fix %s
|
|
|
|
dump 1 %s custom %s ${sname}.%s.lammpstrj id type x y z spin radius
|
|
dump 2 %s custom %s ${sname}.%s.xyz
|
|
|
|
run %s
|
|
|
|
unfix 1
|
|
undump 1
|
|
undump 2
|
|
"""
|
|
|
|
task={'single_pt':single_pt,'minimize':minimize,'dynamics':dynamics}
|
|
|
|
q2m={1:'1.007940',2:'4.002602',3:'6.941000',4:'9.012182',5:'10.811000',6:'12.010700',7:'14.006700',8:'15.999400',
|
|
9:'18.9984032',10:'20.179700',11:'22.98976928',12:'24.305000',13:'26.9815386',14:'28.085500',15:'30.973762',
|
|
16:'32.065000',17:'35.453000',18:'39.948000'}
|
|
|
|
def generate_lammps_input(infile):
|
|
|
|
# Defaults values
|
|
ensemble={"nve":"1 %s nve/eff",'nvt':"1 %s nvt/eff %s %s %s %s",'npt':"1 %s npt/eff %s %s %s %s %s %s"}
|
|
boundary="f f f"
|
|
xbound="-1000.000 1000.0 xlo xhi\n"
|
|
ybound="-1000.000 1000.0 ylo yhi\n"
|
|
zbound="-1000.000 1000.0 zlo zhi\n"
|
|
cutoff=1000.0
|
|
period="1"
|
|
emass=0
|
|
vels=""
|
|
|
|
datafile=open("data."+infile[:-4],'w')
|
|
scriptfile=open("in."+infile[:-4],'w')
|
|
|
|
print "Reading %s ... [WAIT]"%infile,
|
|
fin = open(infile,'r')
|
|
lines = fin.xreadlines()
|
|
print 7*"\b"+"[DONE]"
|
|
|
|
numnuclei=0
|
|
numelec=0
|
|
nuclei={}
|
|
electrons={}
|
|
masses=[]
|
|
massstr="Masses\n\n"
|
|
types=1
|
|
q2type={}
|
|
Tflag=False # Default ensemble is NVE
|
|
steps='1000'
|
|
|
|
print "Extracting run parameters from %s ... "%(infile),
|
|
for line in lines:
|
|
|
|
# 1st level keywords
|
|
if line.find("@params")==0:
|
|
flag='params'
|
|
continue
|
|
elif line.find("@nuclei")==0:
|
|
flag='nuclei'
|
|
continue
|
|
elif line.find("@electrons")==0:
|
|
flag='electrons'
|
|
continue
|
|
elif line.find("@nuc_velocities")==0:
|
|
flag='n_vels'
|
|
continue
|
|
elif line.find("@elec_velocities")==0:
|
|
flag='e_vels'
|
|
continue
|
|
elif line.find("@nuc_masses")==0:
|
|
flag='n_mass'
|
|
continue
|
|
elif line.find("@elec_masses")==0:
|
|
flag='e_mass'
|
|
continue
|
|
elif line.find("@restraints")==0:
|
|
flag='restraints'
|
|
continue
|
|
|
|
# 2nd level keywords
|
|
if flag=='params':
|
|
if line.find("calc")>=0:
|
|
op=line.split()[2]
|
|
if line.find("print_every")>=0:
|
|
period=line.split()[2]
|
|
if line.find("num_steps")>=0:
|
|
steps=line.split()[2]
|
|
if line.find("min_freeze")>=0:
|
|
setforce="velocity\t% set 0.0 0.0 0.0\nfix\tfreeze %s setforce 0.0 0.0 0.0"%(line.split()[2],line.split()[2])
|
|
if line.find("thermostat")>=0:
|
|
tflag=True
|
|
#ensemble="fix\t1 all nvt/eff "
|
|
if line.find("start_temperature")>=0:
|
|
Tstart=line.split()[2]
|
|
#ensemble+=Tstart
|
|
if line.find("end_temperature")>=0:
|
|
Tstop=line.split()[2]
|
|
#ensemble+=Tstop
|
|
if line.find("andersen_coupling")>=0 or line.find("nose_hoover_coupling")>=0:
|
|
Tdamp=line.split()[2]
|
|
#ensemble+=Tdamp
|
|
if line.find("dt")>=0:
|
|
dt=line.split()[2]
|
|
if line.find("electron_mass")>=0:
|
|
emass=line.split()[2]
|
|
if line.find("adaptive_step_size")>=0:
|
|
continue
|
|
if line.find("adaptive_energy")>=0:
|
|
continue
|
|
if line.find("e_field_freq")>=0:
|
|
continue
|
|
if line.find("e_field_packet_duration")>=0:
|
|
continue
|
|
if line.find("e_field")>=0:
|
|
field=line.split()[2:5]
|
|
efield="fix\field all efield %s %s %s"%(field[0],field[1],field[2])
|
|
if line.find("e_field_packet_duration")>=0:
|
|
continue
|
|
if line.find("set_limit")>=0:
|
|
continue # need to add this contraint
|
|
if line.find("set_limit_stiffness")>=0:
|
|
continue
|
|
if line.find("output_position")>=0:
|
|
dump_pos="dump\t1 all custom %s ${sname}.lammpstrj id type x y z spin radius "%(period)
|
|
if line.find("output_velocities")>=0:
|
|
dump_pos+="vx vy vz "
|
|
if line.find("output_energy_forces")>=0:
|
|
dump_pos="compute\tenergy all pe/atom\n"+dump_pos
|
|
dump_pos+="c_energy fx fy fz\n"
|
|
if line.find("output_restart")>=0:
|
|
restart="restart\t%s ${sname}.restart1 ${sname}.restart2"%(period)
|
|
if line.find("output_restraints")>=0:
|
|
continue
|
|
if line.find("ewald_re_cutoff")>=0 or line.find("ewald_autoset")>=0 or line.find("ewald_log_precision")>=0 or line.find("ewald_max_re")>=0 or \
|
|
line.find("ewald_r_cutoff")>=0 or line.find("ewald_k_cutoff")>=0 or line.find("ewald_nuc_r")>=0:
|
|
continue
|
|
if line.find("periodic")>=0:
|
|
bounds=line.split()[2]
|
|
if bounds=="True": boundary="p p p"
|
|
elif bounds=="minimage_x": boundary="p f f"
|
|
elif bounds=="minimage_xy": boundary="p p f"
|
|
elif bounds=="minimage_y": boundary="f p f"
|
|
elif bounds=="minimage_xyz": boundary="p p p"
|
|
elif bounds=="minimage_z": boundary="f f p"
|
|
if line.find("x_bound")>=0:
|
|
xbnds=line.split()[2:4]
|
|
xbound="%s %s xlo xhi\n"%(xbnds[0],xbnds[1])
|
|
if line.find("y_bound")>=0:
|
|
ybnds=line.split()[2:4]
|
|
ybound="%s %s ylo yhi\n"%(ybnds[0],ybnds[1])
|
|
if line.find("z_bound")>=0:
|
|
zbnds=line.split()[2:4]
|
|
zbound="%s %s zlo zhi\n"%(zbnds[0],zbnds[1])
|
|
if line.find("taper_cutoff")>=0:
|
|
cutoff=line.split()[2]
|
|
continue
|
|
|
|
if flag=='nuclei' and len(line)>1:
|
|
numnuclei+=1
|
|
ln=line.split()
|
|
np=' '.join(ln[0:3])
|
|
q=ln[3]
|
|
m=q2m[int(float(q))]
|
|
if m not in masses:
|
|
masses.append(m)
|
|
massstr+="%d %s\n"%(types,m)
|
|
q2type[q]=types
|
|
types+=1
|
|
nuclei[numnuclei]=[np,q]
|
|
continue
|
|
|
|
if flag=='electrons' and len(line)>1:
|
|
numelec+=1
|
|
ln=line.split()
|
|
ep=' '.join(ln[0:3])
|
|
spin=ln[3]
|
|
radius=ln[4]
|
|
electrons[numelec]=[ep,spin,radius]
|
|
if numelec==1:
|
|
if emass!=0: massstr+="%d %s\n\n"%(types,emass) # electron mass=1
|
|
else: massstr+="%d 1.000000\n\n"%(types)
|
|
continue
|
|
|
|
if flag=='n_vels' and len(line)>1:
|
|
vels+=line+" 0.0"
|
|
continue
|
|
|
|
if flag=='e_vels' and len(line)>1:
|
|
ln=line.split()
|
|
ln[0]=ln[0]+numnuclei
|
|
vels+=ln[0]+" "+ln[1]+" "+ln[2]+" "+ln[3]+" "+ln[4]+"\n"
|
|
continue
|
|
|
|
if flag=='n_mass' and len(line)>1:
|
|
print "Setting nuclear masses is unsupported\n"
|
|
continue
|
|
|
|
if flag=='e_mass' and len(line)>1:
|
|
print "Setting electron masses is unsupported\n"
|
|
continue
|
|
|
|
print "\bDone"
|
|
|
|
# Build data file
|
|
print "Writing datafile to %s ... "%('data.'+infile),
|
|
sys.stdout.flush()
|
|
print "\b"*19+"General section ",
|
|
datafile.writelines("Created using cfg2lammps (c) AJB-2009\n\n%d atoms\n%d atom types\n\n%s%s%s\n"%(numnuclei+numelec,types,xbound,ybound,zbound))
|
|
print "\b"*19+"Masses section ",
|
|
datafile.writelines(massstr)
|
|
print "\b"*19+"Atoms section ",
|
|
datafile.writelines("Atoms\n\n")
|
|
for n in range(numnuclei):
|
|
datafile.writelines("%d %d %s %s 0 0.0\n"%(n+1,q2type[nuclei[n+1][1]],nuclei[n+1][0],nuclei[n+1][1]))
|
|
for e in range(numelec):
|
|
datafile.write("%d %d %s 0.0 %s %s\n"%(e+numnuclei+1,types,electrons[e+1][0],electrons[e+1][1],electrons[e+1][2]))
|
|
print "\b"*19+"Velocities section\n",
|
|
datafile.writelines(vels)
|
|
datafile.writelines("\n")
|
|
print "DONE .... GOODBYE !!"
|
|
datafile.close()
|
|
|
|
# Build input script
|
|
import datetime
|
|
scriptfile.writelines(general%(datetime.date.today(),infile[:-4],boundary,cutoff,period))
|
|
if op=='minimize':
|
|
scriptfile.writelines(minimize%('all',period,'all',period,steps,'10000'))
|
|
#%(group,period,group,period,iterations,fcalls)
|
|
elif op=='single_pt':
|
|
scriptfile.writelines(single_pt%())
|
|
elif op=='dynamics':
|
|
if Tflag==True:
|
|
scriptfile.writelines(dynamics%('NVT',dt,ensemble['nvt']%('all',Tstart,Tstop,Tdamp,''),'all',period,'nvt','all',period,'nve',steps))
|
|
#%(ensemble,dt,group,ensemble%(group,tstart,tstop,tdamp,options))
|
|
else:
|
|
scriptfile.writelines(dynamics%('NVE',dt,ensemble['nve']%('all'),'all',period,'nve','all',period,'nve',steps))
|
|
#%(ensemble,dt,group,ensemble%(group))
|
|
scriptfile.writelines("\n")
|
|
|
|
if __name__ == '__main__':
|
|
|
|
# set defaults
|
|
|
|
# check for input:
|
|
opts, argv = getopt(sys.argv[1:], 'h')
|
|
|
|
# if no input, print help and exit
|
|
if len(argv) != 1:
|
|
printHelp()
|
|
sys.exit(1)
|
|
else:
|
|
infile=argv[0]
|
|
|
|
# read options
|
|
for opt, arg in opts:
|
|
if opt == '-h': # -h: print help
|
|
printHelp()
|
|
|
|
generate_lammps_input(infile)
|