git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@5712 f3b2605a-c512-4ea7-a41b-209d697bcdaa

tools/eff/README

@@ -5,9 +5,12 @@ generating LAMMPS structure files:
 - h2.pl: rectangular lattice of hydrogen atoms
 - Diamond.pl: diamond A4 box
 - Li-hydride: Lithium hydride solid box
-- Li-solid: Lithium solid box
+- Li-solid: Lithium solid box (Li-solid-angstromg produces data file in Angstroms)
 - Uniform-electron-gas.pl: uniform electron gas on an NaCl lattice
 
+- bohr2ang.py: converts an eFF structure file described in bohr to 
+  its corresponding version in Angstroms
+
 And other useful scripts for processing pEFF related information are
 included, such as:
 

tools/eff/VMD-input.py

@@ -12,8 +12,8 @@ ajaramil@wag.caltech.edu
 Project: pEFF
 Version: August 2009
 
-Usage: python VMD-input.py lammps_dump_filename radii_column_number
+Usage: python VMD-input.py lammps_dump_filename radii_column xpos_column
-Example: python VMD-input.py dump.lammpstrj 6
+Example: python VMD-input.py dump.lammpstrj 5 6
 
 1. Extracts the electron radii from a lammps trajectory dump into %s.out
 2. Creates %s.xyz file
@@ -28,20 +28,31 @@ def printHelp(input):
   Info%(input,input,input)
 
 if __name__ == '__main__':
-
+    
    # if no input, print help and exit
     if len(sys.argv) < 2:
-        print "Usage: python VMD-input.py lammps_dump_filename radii_column_number\n"
+        print "Usage: python VMD-input.py lammps_dump_filename radii_column xpos_column\n"
         sys.exit(1)
     else:
         infile=sys.argv[1]
 
+    workdir=os.getcwd()
+    tools=sys.argv[0].split('VMD-input.py')[0]
     # set defaults
     outfile = infile.split('.')[0]
-    if len(sys.argv) == 2:
+    print sys.argv
+    if len(sys.argv) == 4:
       column = int(sys.argv[2])
-    else:
+      xpos = int(sys.argv[3])
-      column=6    # default = radius for dump -> id type x y z spin radius
+      print "Assuming xpos=eradius+1"
+    elif len(sys.argv) == 3:
+      column = int(sys.argv[2])
+      xpos = column+1
+      print "Assuming xpos=eradius+1"
+    elif len(sys.argv) == 2:
+      column=5    # default = radius for dump -> id type q spin eradius x y z
+      xpos=6
+    else: print "Incorrect number of arguments"
 
     # check for input:
     opts, argv = getopt(sys.argv[1:], 'c:o:ha')
@@ -54,12 +65,12 @@ if __name__ == '__main__':
           outfile=arg
         if opt == '-c':         # select column from lammpstrj file to tabulate
           column=int(arg)
-
+    print column,xpos
     makeradii(infile,outfile+".out",column,True)
-    lmp2xyz(infile,outfile+".xyz")
+    lmp2xyz(infile,outfile+".xyz",xpos)
     print "Creating %s file ..."%(outfile+".vmd")
-    os.system("cat %s | sed 's/xyzfile/%s/' > %s"%("radii.vmd",outfile+".xyz","temp"))
+    os.system("cat %s | sed 's/xyzfile/%s/' > %s"%(tools+"radii.vmd",outfile+".xyz","temp"))
-    os.system("cat %s | sed 's/radiifile/%s/' > %s; rm temp"%("temp",outfile+".out",outfile+".vmd"))
+    os.system("cat %s | sed 's/radiifile/%s/' > %s; rm temp"%("temp",outfile+".out",workdir+'/'+outfile+".vmd"))
     print "Done !! (you can now source %s using VMD's console) \n"%(outfile+".vmd")
 
     print "NOTE: In VMD, set graphics representation for electrons to transparency,"

tools/eff/cfg2lammps.py

@@ -12,6 +12,8 @@ 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
+
+12/2010: Added support for fixed-core and pseudo-core structures
 """
 
 # import essentials: 
@@ -38,15 +40,31 @@ units		electron
 newton		on
 boundary	%s
 
-atom_style      hybrid charge electron
+atom_style      electron
 
 read_data       data.${sname}
 
 pair_style      eff/cut %s
 pair_coeff      * *
 
+compute         energies all pair eff/cut
+variable        eke equal c_energies[1]
+variable        epauli equal c_energies[2]
+variable        estatics equal c_energies[3]
+variable        errestrain equal c_energies[4]
+
+communicate     single vel yes
+
+compute         peratom all stress/atom
+compute         p all reduce sum c_peratom[1] c_peratom[2] c_peratom[3]
+variable        press equal -(c_p[1]+c_p[2]+c_p[3])/(3*vol)
+
+compute         effTemp all temp/eff
+compute         effPress all pressure effTemp
+
 thermo          %s
-thermo_style    multi
+thermo_style    custom step etotal pe ke v_eke v_epauli v_estatics v_errestrain temp press v_press 
+thermo_modify   temp effTemp press effPress
 """
 #%(date,name,boundary,cutoff,period)
 
@@ -55,7 +73,7 @@ minimize="""
 
 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
+dump            2 %s custom %s ${sname}.min.lammpstrj id type q spin eradius x y z fx fy fz erforce
 min_modify      line quadratic
 minimize        0 1.0e-5 %s %s
 
@@ -77,7 +95,7 @@ timestep        %s
 
 fix             %s
 
-dump            1 %s custom %s ${sname}.%s.lammpstrj id type x y z spin radius
+dump            1 %s custom %s ${sname}.%s.lammpstrj id type q spin eradius x y z
 dump		2 %s custom %s ${sname}.%s.xyz
 
 run             %s
@@ -114,8 +132,10 @@ def generate_lammps_input(infile):
     lines = fin.xreadlines()
     print 7*"\b"+"[DONE]"
     
+    numcores=0
     numnuclei=0
     numelec=0
+    cores={}
     nuclei={}
     electrons={}
     masses=[]
@@ -132,6 +152,9 @@ def generate_lammps_input(infile):
       if line.find("@params")==0:
         flag='params'
         continue
+      elif line.find("@cores")==0:
+        flag='cores'
+        continue
       elif line.find("@nuclei")==0:
         flag='nuclei'
         continue
@@ -198,7 +221,7 @@ def generate_lammps_input(infile):
         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)
+          dump_pos="dump\t1 all custom %s ${sname}.lammpstrj id type q spin eradius x y z "%(period)
         if line.find("output_velocities")>=0:
           dump_pos+="vx vy vz "
         if line.find("output_energy_forces")>=0:
@@ -232,6 +255,22 @@ def generate_lammps_input(infile):
           cutoff=line.split()[2]            
         continue
 
+      if flag=='cores' and len(line)>1:
+        numcores+=1
+        ln=line.split()
+        nc=' '.join(ln[0:3])
+        q=ln[3]
+        spin='3'
+        radius=ln[4]
+        m=q2m[int(float(q))]
+        if m not in masses:
+          masses.append(m)
+          massstr+="%d %s\n"%(types,m)
+          q2type[q]=types
+          types+=1
+        cores[numcores]=[nc,q,spin,radius]
+        continue
+
       if flag=='nuclei' and len(line)>1:
         numnuclei+=1
         ln=line.split()
@@ -282,15 +321,17 @@ def generate_lammps_input(infile):
     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))
+    datafile.writelines("Created using cfg2lammps (c) AJB-2009\n\n%d atoms\n%d atom types\n\n%s%s%s\n"%(numcores+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(numcores):
+      datafile.writelines("%d %d %2.2f %s %s %s\n"%(n+1,q2type[cores[n+1][1]],float(cores[n+1][1]),cores[n+1][2],cores[n+1][3],cores[n+1][0]))
     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]))
+      datafile.writelines("%d %d %2.2f 0 0.0 %s\n"%(n+numcores+1,q2type[nuclei[n+1][1]],float(nuclei[n+1][1]),nuclei[n+1][0]))
     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]))
+      datafile.write("%d %d 0.0 %s %s %s\n"%(e+numnuclei+numcores+1,types,electrons[e+1][1],electrons[e+1][2],electrons[e+1][0]))
     print "\b"*19+"Velocities section\n",
     datafile.writelines(vels)
     datafile.writelines("\n")

tools/eff/lmp2radii.py

@@ -11,9 +11,9 @@ Version: August 2009
 
 Extracts the electron radii from a lammps trajectory dump of style custom:
 
-dump    1 all custom period dump_file id type x y z spin radius ...
+dump    1 all custom period dump_file id type q spin eradius x y z...
 
-NOTE: The radius must be the 6th column per trajectory entry in the dump file
+NOTE: The radius must be the 5th column per trajectory entry in the dump file
 
 """
 
@@ -45,7 +45,7 @@ def makeradii(infile):
     tmp=open("atoms",'r')
     atoms=int(tmp.readlines()[1].split()[0])
     tmp.close()
-    os.system("rm -rf frames atoms")
+    os.system("rm -rf frames atoms lines")
     arry=numpy.zeros((atoms,frames),dtype=float)
     framecnt=0
     header=9
@@ -59,8 +59,8 @@ def makeradii(infile):
       elif (i >= lo) and (i < hi):
         lparse=line.split()
         id=int(lparse[0])
-        r=float(lparse[6])
+        r=float(lparse[4])
-        if (r!=0): 
+        if (r!=0.0): 
           arry[id-1][framecnt]=r
           if (framecnt==0): ecount+=1
         if (i==lo+1): 

tools/eff/lmp2radii.pyx

@@ -11,9 +11,9 @@ Version: August 2009
 
 Extracts the electron radii from a lammps trajectory dump of style custom:
 
-dump    1 all custom period dump_file id type x y z spin radius ...
+dump    1 all custom period dump_file id type q spin eradius x y z...
 
-NOTE: The radius must be the 6th column per trajectory entry in the dump file
+NOTE: The radius must be the 5th column per trajectory entry in the dump file
 
 """
 
@@ -59,7 +59,7 @@ def makeradii(infile):
       elif (i >= lo) and (i < hi):
         lparse=line.split()
         id=int(lparse[0])
-        r=float(lparse[6])
+        r=float(lparse[4])
         if (r!=0): 
           arry[id-1][framecnt]=r
           if (framecnt==0): ecount+=1
@@ -87,7 +87,7 @@ def makeradii(infile):
           fout.writelines("%f\t"%(arry[a][f]))
         fout.writelines("\n")
     print
-    print "DONE .... GOODBYE !!"
+    print "Done !! (generated radii/frame table) \n"
     fout.close()
     fin.close()
 

tools/eff/lmp2radii_col.py

@@ -11,7 +11,7 @@ Version: August 2009
 
 Extracts the electron radii from a lammps trajectory dump of style custom:
 
-dump    1 all custom period dump_file id type x y z spin radius ...
+dump    1 all custom period dump_file id type spin eradius x y z ...
 
 NOTE: The radius must be the "column" per trajectory entry in the dump file
 
@@ -59,7 +59,7 @@ def makeradii(infile,outfile,column,True):
       elif (i >= lo) and (i < hi):
         lparse=line.split()
         id=int(lparse[0])
-        r=float(lparse[column])
+        r=float(lparse[column-1])
         if (r!=0): 
           arry[id-1][framecnt]=r
           if (framecnt==0): ecount+=1
@@ -96,7 +96,7 @@ if __name__ == '__main__':
     # set defaults
     outfile = ""
     flag_all = False
-    column=6    # default = radius
+    column=5    # default = radius
 
     # check for input:
     opts, argv = getopt(sys.argv[1:], 'c:o:ha')

tools/eff/lmp2xyz.py

@@ -15,15 +15,35 @@ Usage: python lmp2xyz.py lammps_dump_filename xyz_filename
 """
 
 import os, sys
-from math import log10
+from math import log10,floor
 from numpy import zeros
 
-mass={"1.00794":"H","4.002602":"He","6.941":"Li","9.012182":"Be","10.811":"B","12.0107":"C","1.00":"Au","0.0005486":"Au"}
+masses={"1.00794":"H","4.002602":"He","6.941":"Li","9.012182":"Be","10.811":"B","12.0107":"C","1.00":"Au","0.0005486":"Au"}
+mass_floor={1:"H",4:"He",6:"Li",9:"Be",10:"B",12:"C",0:"Au",28:"Si"}
 
-def lmp2xyz(lammps,xyz):
+def lmp2xyz(lammps,xyz,xpos):
   print "\nGenerating %s file"%(xyz)
   fin=open(lammps,'r')
   fout=open(xyz,'w')
+  data=raw_input("Do you have a corresponding data file? please enter filename or 'n': ")
+  count=1
+  if data!='n': 
+    dataf=open(data,'r')
+    datafile=dataf.readlines()
+    dataf.close()
+    for line in datafile:
+      if line.find("atom types")>=0:
+        numtypes=int(line.split()[0])
+        mass=zeros(numtypes,dtype=float)
+      elif line.find("Masses")>=0:
+        count+=1+datafile.index(line)
+      elif line.find("Atoms")>=0:
+        break
+    for i in range(numtypes):
+      mass[i]=float(datafile[count].split()[1])
+      count+=1
+  else: 
+    print "\nWill continue without a data file specification"
   header=9
   lines=fin.readlines()
   numatoms=lines[3].split()[0]
@@ -31,13 +51,15 @@ def lmp2xyz(lammps,xyz):
   tmp=open('lines','r')
   tlines=tmp.readline()
   tmp.close()
+  os.system("rm lines")
   flines=int(tlines.split()[0])
   snaps=flines/(int(numatoms)+header)
   countsnap=1
-  coords=zeros((int(numatoms),4),dtype=float)
+  if data!='n': coords={}
+  else: coords=zeros((int(numatoms),4),dtype=float)
   sys.stdout.write("Writing [%d]: "%(snaps))
   sys.stdout.flush()
-  read_atoms=1
+  read_atoms=0
   for line in lines:
     if line.find('ITEM: TIMESTEP')==0:
       read_atom_flag=False
@@ -53,14 +75,23 @@ def lmp2xyz(lammps,xyz):
       read_atoms+=1
       parse=line.split()
       if parse[0]!="":
-        coords[int(parse[0])-1][0]=int(parse[1])
+#        print [mass_floor[int(floor(mass[int(parse[1])-1]))],float(parse[xpos-1]),float(parse[xpos]),float(parse[xpos+1])]
-        coords[int(parse[0])-1][1]=float(parse[2])
+        if data!='n': 
-        coords[int(parse[0])-1][2]=float(parse[3])
+          if mass[int(parse[1])-1]==1.0:
-        coords[int(parse[0])-1][3]=float(parse[4])
+            type='Au'
+          else:
+            type=mass_floor[int(floor(mass[int(parse[1])-1]))]
+          coords[int(parse[0])-1]=[type,float(parse[xpos-1]),float(parse[xpos]),float(parse[xpos+1])]
+        else: 
+          coords[int(parse[0])-1][0]=int(parse[1])
+          coords[int(parse[0])-1][1]=float(parse[xpos-1])
+          coords[int(parse[0])-1][2]=float(parse[xpos])
+          coords[int(parse[0])-1][3]=float(parse[xpos+1])
       if read_atoms==int(numatoms):
         read_atoms=0
         for i in range(int(numatoms)):
-          fout.writelines("%d %2.4f %2.4f %2.4f\n"%(coords[i][0],coords[i][1],coords[i][2],coords[i][3]))
+          if data!='n': fout.writelines("%s %2.4f %2.4f %2.4f\n"%(coords[i][0],coords[i][1],coords[i][2],coords[i][3]))
+          else: fout.writelines("%d %2.4f %2.4f %2.4f\n"%(coords[i][0],coords[i][1],coords[i][2],coords[i][3]))
 
   print "\nDone converting to xyz!!\n"
   fin.close()
@@ -76,6 +107,8 @@ if __name__ == '__main__':
 
     inputfile=sys.argv[1]
     outfile=sys.argv[2]
-
+    if len(sys.argv)==4:
-    lmp2xyz(inputfile,outfile.split()[0])
+      xpos=sys.arv[3]-1
+    else: xpos=5
+    lmp2xyz(inputfile,outfile.split()[0],xpos)
 

tools/eff/radii.vmd

@@ -171,7 +171,7 @@ proc do_radii {args} {
   foreach elec $data_array(0) r $data_array($fr) {
     set s [atomselect $molid "index [expr {$elec -1}]"]
     #set nr [expr {exp($r)}]
-    set nr $r
+    set nr [expr $r*1.5]
     $s set radius $nr
     $s delete
     #puts stderr "$fr $elec $nr"
@@ -192,9 +192,6 @@ global data_array molid
 set xyz xyzfile
 set data radiifile
 
-# switch default rep to VDW.
-mol default style VDW
-
 # load nuclear and electron xyz trajectory 
 #set xyz [lindex $::argv 0]
 
@@ -202,7 +199,35 @@ mol default style VDW
 #set datafile [lindex $::argv 1]
 
 set molid [mol new $xyz waitfor all]
-mol modstyle 0 [molinfo top] VDW 1.0 32.0
+#mol modstyle 0 [molinfo top] VDW 0.2 32.0
+
+set s [atomselect $molid "all"]
+$s set radius 0.1
+$s delete
+
+color Display Background white
+display projection orthographic
+
+# switch default rep to VDW.
+mol default style VDW
+mol delrep 0 $molid
+mol selection "name H He Li Be B C Si"
+mol addrep $molid
+mol modstyle 0 $molid VDW 2.0 32.0
+# spin +1
+mol material Transparent
+mol selection "name Au"
+mol addrep $molid
+mol modstyle 1 $molid VDW 0.4 32.0
+# spin -1
+#mol material Transparent
+mol selection "name Tl"
+mol addrep $molid
+mol modstyle 2 $molid VDW 0.4 32.0
+mol material Opaque
+mol selection "name H He Li Be B C Si"
+mol addrep $molid
+mol modstyle 3 $molid DynamicBonds 1.6 0.1 32
 
 puts "Starting ..."
 readFile $data