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

tools/eff/README

+The directory structure_generators includes several scripts for
+generating LAMMPS structure files:
+
+- Be-solid.pl: beryllium solid box
+- h2.pl: rectangular lattice of hydrogen atoms
+- Diamond.pl: diamond A4 box
+- Li-hydride: Lithium hydride solid box
+- Li-solid: Lithium solid box
+- Uniform-electron-gas.pl: uniform electron gas on an NaCl lattice
+
+And other useful scripts for processing pEFF related information are
+included, such as:
+
+- cfg2lammps.py: Python script for converting an eff cfg file into a
+  LAMMPS data/script file pair
+
+- lmp2radii.py/pyx: Python/Cython scripts for post-processing a lammps
+  trajectory to extract electron radii/frame.
+
+Note: the corresponding .c, and .so files are the c source and binary
+library (loadable Python module) created by Cython automatically when
+compiling lmp2radii.pyx in place with Cython (python setup.py
+build_ext --inplace).
+
+- radii.vmd: a TCL script for adding radial changes per trajectory
+  frame to an xyz LAMMPS trajectory of an pEFF run.
+
+- lmp2any.py: Python scipt for extracting quantities from a custom
+  lammps dump
+
+- lmp2radii-col.py: Same as lmp2radii.py, but takes a column
+  descriptor for the radius
+
+- VMD-input.py: Automatically calls the necessary scripts to produce a
+  VMD ready script that loads variable radii into VMD
+
+NOTE: you must set the graphical representation in VMD to represent
+electrons using transparency (this requires selecting and applying the
+corresponding atom types).
+
+For further details see the descriptors in each file, or contact:
+Andres Jaramillo-Botero: ajaramil@wag.caltech.edu
+
+Acknowledgments: Thanks to Axel Kohlmeyer (Temple Univ) for his help
+with VMD.

tools/eff/VMD-input.py

+#!/usr/local/bin/python-2.5/bin/python
+
+import sys, os
+from getopt import gnu_getopt as getopt
+
+Info="""
+Module name: VMD-input.py 
+
+Author: (c) Andres Jaramillo-Botero
+California Institute of Technology
+ajaramil@wag.caltech.edu
+Project: pEFF
+Version: August 2009
+
+Usage: python VMD-input.py lammps_dump_filename radii_column_number
+Example: python VMD-input.py dump.lammpstrj 6
+
+1. Extracts the electron radii from a lammps trajectory dump into %s.out
+2. Creates %s.xyz file
+3. Creates appropriate %.vmd file which can be sourced using TCL/TK in VMD
+
+"""
+
+from lmp2radii_col import makeradii 
+from lmp2xyz import lmp2xyz
+ 
+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"
+        sys.exit(1)
+    else:
+        infile=sys.argv[1]
+
+    # set defaults
+    outfile = infile.split('.')[0]
+    if len(sys.argv) == 2:
+      column = int(sys.argv[2])
+    else:
+      column=6    # default = radius for dump -> id type x y z spin radius
+
+    # check for input:
+    opts, argv = getopt(sys.argv[1:], 'c:o:ha')
+
+    # read options
+    for opt, arg in opts:
+        if opt == '-h':             # -h: print help
+          Info%(input,input,input)
+        if opt == '-o':         # output file name
+          outfile=arg
+        if opt == '-c':         # select column from lammpstrj file to tabulate
+          column=int(arg)
+
+    makeradii(infile,outfile+".out",column,True)
+    lmp2xyz(infile,outfile+".xyz")
+    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/radiifile/%s/' > %s; rm temp"%("temp",outfile+".out",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,"
+    print "and change the atom types in the xyz file according to your values,"
+    print "for simplicity, they are set using the same mass sequence definition\nfrom your lammps data file\n"

tools/eff/cfg2lammps.py

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

tools/eff/h2.vmd

+# module radii.vmd
+# December, 2009 -(c)- Andres Jaramillo-Botero
+# Script to load variable changing radii onto a pEFF lammpstrj file
+
+# radii.vmd --
+#    Script to read and change electron radii in vmd dynamics traj
+#
+ 
+# openFile --
+#    Open the file and start looking for data
+#
+# Arguments:
+#    filename    Name of the file to read
+#
+# Result:
+#    infile      Handle to the opened file
+#
+# Side effects:
+#    A file in question is opened. Lines containing a * or# as
+#    the first non-blank character are considered comments. The
+#    first line without this is considered to be a line with the
+#    names of the columns.
+#
+proc openFile {filename} {
+    set infile [open $filename "r"]
+    return $infile
+}
+
+# readNames --
+#    Read the names of the columns
+#
+# Arguments:
+#    infile      Handle to the file
+#
+# Result:
+#    names       List of names, the number indicates the number of
+#                the snapshot frame
+#
+proc readNames {infile} {
+  #
+  # Skip the header - if any
+  #
+    set pos    0
+    while { [gets $infile line] >= 0 } {
+       if { [regexp {[ \t]*[*#]} $line] } {
+          incr pos
+       } else {
+          break
+       }
+    }
+    seek $infile 0 start
+    while { $pos > 0 } {
+       gets $infile line
+       incr pos -1
+    }
+
+  #
+  # Read the line with the column names
+  #
+    gets $infile line
+
+  # Force the line to be interpreted as a list
+    set nocols [llength $line]
+
+    return $line
+}
+
+# readData --
+#    Read the data per line
+#
+# Arguments:
+#    infile      Handle to the file
+#
+# Result:
+#    values      List of values, representing each column.
+#                A list of length zero indicates the end of file
+#
+proc readData {infile} {
+    while { [gets $infile values] == 0 } { ;# Just go on - skip empty lines }
+
+    set nocols [llength $values]
+
+    return $values
+}
+
+# readFile --
+#    Read the file and store the data in a (global) array
+#
+# Arguments:
+#    filename    Name of the file
+#
+# Result:
+#    None
+#
+# Side effects:
+#    Filled array, ready for display
+#
+proc readFile {filename} {
+    global data_array
+
+    set infile [openFile $filename]
+
+    set data_array(names) [readNames $infile]
+    set i 0
+    foreach name $data_array(names) {
+       set data_array($i) {}
+       incr i
+    }
+
+    while 1 {
+       set values [readData $infile]
+       if { [llength $values] > 0 } {
+          set i 0
+          foreach value $values {
+             lappend data_array($i) $value
+             incr i
+          }
+       } else {
+          break
+       }
+    }
+}
+
+# makeXYData --
+#    Make a list useable by frame-electron
+#
+# Arguments:
+#    xindex     Index of frame data
+#    yindex     Index of electron data
+#
+# Result:
+#    None
+#
+# Side effects:
+#    A dataset for changing the electron radii, per trajectory frame
+#
+proc makeXYData {xindex yindex} {
+    global data_array
+
+    set xydata {}
+    foreach x $data_array($xindex) y $data_array($yindex) {
+        lappend xydata $x $y
+    }
+    return $xydata
+}
+
+proc returnXYPair {x y} {
+   global data_array
+
+   return [list $data_array($x) $data_array($y)]
+}
+
+# do_radii --
+#    Changes the radii of electrons per trajectory frame
+#
+# Arguments:
+#    xindex     Index of frame data
+#    yindex     Index of electron data
+#
+# Result:
+#    prints the frame:atomID:radius
+#
+proc do_radii {args} {
+  global molid data_array
+
+  #set n [molinfo $molid get numatoms]
+  set f [molinfo $molid get frame]
+ 
+  set fr [expr {$f+1}]
+
+  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
+    $s set radius $nr
+    $s delete
+    #puts stderr "$fr $elec $nr"
+  }
+}
+
+# main --
+#    Main control flow
+#
+
+# Check input arguments
+#for {set i 0} {$i<[llength $argv]} {incr i}
+#  puts " - $i: [lindex $argv $i]"
+
+global data_array molid
+
+# Set input files manually
+set xyz h2.xyz
+set data h2.out
+
+# switch default rep to VDW.
+mol default style VDW
+
+# load nuclear and electron xyz trajectory 
+#set xyz [lindex $::argv 0]
+
+# load electron radii information for trajectory
+#set datafile [lindex $::argv 1]
+
+set molid [mol new $xyz waitfor all]
+mol modstyle 0 [molinfo top] VDW 1.0 32.0
+
+puts "Starting ..."
+readFile $data
+puts "Read datafile and created array of radii ..."
+
+puts "Visualize trajectory" 
+trace variable vmd_frame($molid) w do_radii
+animate goto start
+do_radii
+  

tools/eff/lmp2data.py

+#!/usr/local/bin/python-2.5/bin/python
+
+Info="""
+Module name: lmp2data.py 
+
+Author: (c) Andres Jaramillo-Botero
+California Institute of Technology
+ajaramil@wag.caltech.edu
+Project: pEFF
+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 ...
+
+NOTE: The radius must be the i'th column per trajectory entry in the dump file
+
+"""
+
+# 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 lmp2data.py test.lammpstrj\n"
+    return
+
+def makeradii(infile,outfile,column,flag_all):
+
+    print "Reading %s ... [WAIT]"%infile,
+    fin = open(infile,'r')
+    lines = fin.xreadlines()
+    print 7*"\b"+"[DONE]"
+    frame=0
+    radii=[]
+    # grep the number of frames and atoms/frame
+    os.system("grep TIMESTEP %s | wc -l > frames; grep -m 1 -A 1 ATOMS %s > atoms; grep -m 1 \"ITEM: ATOMS\" %s > params"%(infile,infile,infile))
+    tmp=open("frames",'r')
+    frames=int(tmp.readline().split()[0])
+    tmp.close()
+    tmp=open("atoms",'r')
+    atoms=int(tmp.readlines()[1].split()[0])
+    tmp.close()
+    tmp=open("params",'r')
+    ids=tmp.readline().split()[2:]
+    os.system("rm -rf frames atoms params")
+    arry=numpy.zeros((atoms,frames),dtype=str)
+    framecnt=0
+    header=9
+    ecount=0
+    if flag_all==True: atom_type="nuclei and electron"
+    else: atom_type="electron"
+    print "Extracting %s %s per frame from %s ...  "%(atom_type,ids[column],infile),
+    for i,line in enumerate(lines):
+      lo=(atoms+header)*framecnt+header
+      hi=lo+atoms
+      if (i<lo): 
+        continue
+      elif (i >= lo) and (i < hi):
+        lparse=line.split()
+        id=int(lparse[0])
+#        r=float(lparse[column-1])
+        r=lparse[column]
+#        if (float(r)!=0):
+        arry[id-1][framecnt]=r
+        print arry[id-1][framecnt],r,raw_input()
+        if (float(r)!=0) and (framecnt==0): ecount+=1
+#        else: arry[id-1][framecnt]=r
+        if (i==lo+1): 
+          sys.stdout.write("%d/%d%s"%(framecnt+1,frames,(int(log10(framecnt+1))+3+int(log10(frames)))*"\b"))
+          sys.stdout.flush()
+      if (i == hi+1): 
+        framecnt+=1
+
+    print
+    if outfile=="": 
+      outfile=infile+'.%s'%(ids[column])
+      fout=open(outfile,'w')
+    else: fout=open(outfile,'w')
+    print "Writing %s/frame table to %s ... "%(ids[column],outfile),
+    sys.stdout.flush()
+
+    for i in range(frames):
+      fout.writelines('\tF'+str(i))
+    fout.writelines("\n")
+    e=1
+    for a in range(atoms):
+      if flag_all==True:
+        sys.stdout.write("%d/%d%s"%(a+1,atoms,(int(log10(a+1))+int(log10(atoms))+3)*"\b"))
+        sys.stdout.flush()
+        fout.writelines("%d\t"%(a+1))
+        for f in range(frames):
+          fout.writelines("%s\t"%(arry[a][f]))
+        fout.writelines("\n")
+      else:
+        if arry[a][0] == 0.0:
+          continue
+        else:
+          sys.stdout.write("%d/%d%s"%(e,ecount,(int(log10(e))+int(log10(ecount))+3)*"\b"))
+          sys.stdout.flush()
+          e+=1
+          fout.writelines("%d\t"%(a+1))
+          for f in range(frames):
+            fout.writelines("%s\t"%(arry[a][f]))
+          fout.writelines("\n")
+
+    print
+    print "DONE .... GOODBYE !!"
+    fout.close()
+    fin.close()
+
+if __name__ == '__main__':
+
+    # set defaults
+    outfile = ""
+    flag_all = False
+    column=6	# default = radius
+
+    # check for input:
+    opts, argv = getopt(sys.argv[1:], 'c:o:ha')
+
+    # 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()
+        if opt == '-o':		# output file name
+          outfile=arg
+        if opt == '-a':		# all nuclii+electrons
+          flag_all=True
+        if opt == '-c':		# select column from lammpstrj file to tabulate
+          column=int(arg)
+
+    makeradii(infile,outfile,column,flag_all)

tools/eff/lmp2radii.py

+#!/usr/local/bin/python-2.5/bin/python
+
+Info="""
+Module name: lmp2radii.py 
+
+Author: (c) Andres Jaramillo-Botero
+California Institute of Technology
+ajaramil@wag.caltech.edu
+Project: pEFF
+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 ...
+
+NOTE: The radius must be the 6th column per trajectory entry in the dump file
+
+"""
+
+# 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 lmp2radii.pyx test.lammpstrj\n"
+    return
+
+def makeradii(infile):
+
+    print "Reading %s ... [WAIT]"%infile,
+    fin = open(infile,'r')
+    lines = fin.xreadlines()
+    print 7*"\b"+"[DONE]"
+    frame=0
+    radii=[]
+    # grep the number of frames and atoms/frame
+    os.system("grep TIMESTEP %s | wc -l > frames; grep -m 1 -A 1 ATOMS %s > atoms"%(infile,infile))
+    tmp=open("frames",'r')
+    frames=int(tmp.readline().split()[0])
+    tmp.close()
+    tmp=open("atoms",'r')
+    atoms=int(tmp.readlines()[1].split()[0])
+    tmp.close()
+    os.system("rm -rf frames atoms")
+    arry=numpy.zeros((atoms,frames),dtype=float)
+    framecnt=0
+    header=9
+    ecount=0
+    print "Extracting electron radii per frame from %s ...  "%(infile),
+    for i,line in enumerate(lines):
+      lo=(atoms+header)*framecnt+header
+      hi=lo+atoms
+      if (i<lo): 
+        continue
+      elif (i >= lo) and (i < hi):
+        lparse=line.split()
+        id=int(lparse[0])
+        r=float(lparse[6])
+        if (r!=0): 
+          arry[id-1][framecnt]=r
+          if (framecnt==0): ecount+=1
+        if (i==lo+1): 
+          sys.stdout.write("%d/%d%s"%(framecnt+1,frames,(int(log10(framecnt+1))+3+int(log10(frames)))*"\b"))
+          sys.stdout.flush()
+      if (i == hi+1): 
+        framecnt+=1
+    print
+    print "Writing radii/frame table to %s ... "%(infile+'.out'),
+    sys.stdout.flush()
+    fout=open(infile+'.out','w')
+    for i in range(frames):
+      fout.writelines('\tF'+str(i))
+    fout.writelines("\n")
+    e=1
+    for a in range(atoms):
+      if arry[a][0] == 0.0: continue
+      else:
+        sys.stdout.write("%d/%d%s"%(e,ecount,(int(log10(e))+int(log10(ecount))+3)*"\b"))
+        sys.stdout.flush()
+        e+=1
+        fout.writelines("%d\t"%(a+1))
+        for f in range(frames):
+          fout.writelines("%f\t"%(arry[a][f]))
+        fout.writelines("\n")
+    print
+    print "Done !! (generated radii/frame table) \n"
+    fout.close()
+    fin.close()
+
+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()
+
+    makeradii(infile)

tools/eff/lmp2radii.pyx

+#!/usr/local/bin/python-2.5/bin/python
+
+Info="""
+Module name: lmp2radii.py 
+
+Author: (c) Andres Jaramillo-Botero
+California Institute of Technology
+ajaramil@wag.caltech.edu
+Project: pEFF
+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 ...
+
+NOTE: The radius must be the 6th column per trajectory entry in the dump file
+
+"""
+
+# 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 lmp2radii.pyx test.lammpstrj\n"
+    return
+
+def makeradii(infile):
+
+    print "Reading %s ... [WAIT]"%infile,
+    fin = open(infile,'r')
+    lines = fin.xreadlines()
+    print 7*"\b"+"[DONE]"
+    frame=0
+    radii=[]
+    # grep the number of frames and atoms/frame
+    os.system("grep TIMESTEP %s | wc -l > frames; grep -m 1 -A 1 ATOMS %s > atoms"%(infile,infile))
+    tmp=open("frames",'r')
+    frames=int(tmp.readline().split()[0])
+    tmp.close()
+    tmp=open("atoms",'r')
+    atoms=int(tmp.readlines()[1].split()[0])
+    tmp.close()
+    os.system("rm -rf frames atoms")
+    arry=numpy.zeros((atoms,frames),dtype=float)
+    framecnt=0
+    header=9
+    ecount=0
+    print "Extracting electron radii per frame from %s ...  "%(infile),
+    for i,line in enumerate(lines):
+      lo=(atoms+header)*framecnt+header
+      hi=lo+atoms
+      if (i<lo): 
+        continue
+      elif (i >= lo) and (i < hi):
+        lparse=line.split()
+        id=int(lparse[0])
+        r=float(lparse[6])
+        if (r!=0): 
+          arry[id-1][framecnt]=r
+          if (framecnt==0): ecount+=1
+        if (i==lo+1): 
+          sys.stdout.write("%d/%d%s"%(framecnt+1,frames,(int(log10(framecnt+1))+3+int(log10(frames)))*"\b"))
+          sys.stdout.flush()
+      if (i == hi+1): 
+        framecnt+=1
+    print
+    print "Writing radii/frame table to %s ... "%(infile+'.out'),
+    sys.stdout.flush()
+    fout=open(infile+'.out','w')
+    for i in range(frames):
+      fout.writelines('\tF'+str(i))
+    fout.writelines("\n")
+    e=1
+    for a in range(atoms):
+      if arry[a][0] == 0.0: continue
+      else:
+        sys.stdout.write("%d/%d%s"%(e,ecount,(int(log10(e))+int(log10(ecount))+3)*"\b"))
+        sys.stdout.flush()
+        e+=1
+        fout.writelines("%d\t"%(a+1))
+        for f in range(frames):
+          fout.writelines("%f\t"%(arry[a][f]))
+        fout.writelines("\n")
+    print
+    print "DONE .... GOODBYE !!"
+    fout.close()
+    fin.close()
+
+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()
+
+    makeradii(infile)

tools/eff/lmp2radii_col.py

+#!/usr/local/bin/python-2.5/bin/python
+
+Info="""
+Module name: lmp2radii-column.py 
+
+Author: (c) Andres Jaramillo-Botero
+California Institute of Technology
+ajaramil@wag.caltech.edu
+Project: pEFF
+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 ...
+
+NOTE: The radius must be the "column" per trajectory entry in the dump file
+
+"""
+
+# 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 lmp2radii.pyx test.lammpstrj\n"
+    return
+
+def makeradii(infile,outfile,column,True):
+
+    print "Reading %s ... [WAIT]"%infile,
+    fin = open(infile,'r')
+    lines = fin.xreadlines()
+    print 7*"\b"+"[DONE]"
+    frame=0
+    radii=[]
+    # grep the number of frames and atoms/frame
+    os.system("grep TIMESTEP %s | wc -l > frames; grep -m 1 -A 1 ATOMS %s > atoms"%(infile,infile))
+    tmp=open("frames",'r')
+    frames=int(tmp.readline().split()[0])
+    tmp.close()
+    tmp=open("atoms",'r')
+    atoms=int(tmp.readlines()[1].split()[0])
+    tmp.close()
+    os.system("rm -rf frames atoms")
+    arry=numpy.zeros((atoms,frames),dtype=float)
+    framecnt=0
+    header=9
+    ecount=0
+    print "Extracting electron radii per frame from %s ...  "%(infile),
+    for i,line in enumerate(lines):
+      lo=(atoms+header)*framecnt+header
+      hi=lo+atoms
+      if (i<lo): 
+        continue
+      elif (i >= lo) and (i < hi):
+        lparse=line.split()
+        id=int(lparse[0])
+        r=float(lparse[column])
+        if (r!=0): 
+          arry[id-1][framecnt]=r
+          if (framecnt==0): ecount+=1
+        if (i==lo+1): 
+          sys.stdout.write("%d/%d%s"%(framecnt+1,frames,(int(log10(framecnt+1))+3+int(log10(frames)))*"\b"))
+          sys.stdout.flush()
+      if (i == hi+1): 
+        framecnt+=1
+    print
+    print "Writing radii/frame table to %s ... "%(infile+'.out'),
+    sys.stdout.flush()
+    fout=open(outfile,'w')
+    for i in range(frames):
+      fout.writelines('\tF'+str(i))
+    fout.writelines("\n")
+    e=1
+    for a in range(atoms):
+      if arry[a][0] == 0.0: continue
+      else:
+        sys.stdout.write("%d/%d%s"%(e,ecount,(int(log10(e))+int(log10(ecount))+3)*"\b"))
+        sys.stdout.flush()
+        e+=1
+        fout.writelines("%d\t"%(a+1))
+        for f in range(frames):
+          fout.writelines("%f\t"%(arry[a][f]))
+        fout.writelines("\n")
+    print
+    print "Done !! (generated radii/frame table) \n"
+    fout.close()
+    fin.close()
+
+if __name__ == '__main__':
+
+    # set defaults
+    outfile = ""
+    flag_all = False
+    column=6    # default = radius
+
+    # check for input:
+    opts, argv = getopt(sys.argv[1:], 'c:o:ha')
+
+    # 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()
+        if opt == '-o':         # output file name
+          outfile=arg
+        if opt == '-a':         # all nuclii+electrons
+          flag_all=True
+        if opt == '-c':         # select column from lammpstrj file to tabulate
+          column=int(arg)
+
+    makeradii(infile,outfile,column,flag_all)
+

tools/eff/lmp2xyz.py

+Info="""
+Module name: lmp2xyz.py 
+
+Author: (c) Andres Jaramillo-Botero
+California Institute of Technology
+ajaramil@caltech.edu
+Project: pEFF
+Version: August 2009
+
+Extracts the xyz from a lammps trajectory dump of style custom:
+dump    1 all custom period dump_file id type x y z spin radius ...
+
+Usage: python lmp2xyz.py lammps_dump_filename xyz_filename
+
+"""
+
+import os, sys
+from math import log10
+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"}
+
+def lmp2xyz(lammps,xyz):
+  print "\nGenerating %s file"%(xyz)
+  fin=open(lammps,'r')
+  fout=open(xyz,'w')
+  header=9
+  lines=fin.readlines()
+  numatoms=lines[3].split()[0]
+  fsize=os.system("wc -l %s> lines"%(lammps))
+  tmp=open('lines','r')
+  tlines=tmp.readline()
+  tmp.close()
+  flines=int(tlines.split()[0])
+  snaps=flines/(int(numatoms)+header)
+  countsnap=1
+  coords=zeros((int(numatoms),4),dtype=float)
+  sys.stdout.write("Writing [%d]: "%(snaps))
+  sys.stdout.flush()
+  read_atoms=1
+  for line in lines:
+    if line.find('ITEM: TIMESTEP')==0:
+      read_atom_flag=False
+      sys.stdout.write("%d "%(countsnap))
+      sys.stdout.flush()
+      fout.writelines("%s\nAtoms\n"%(numatoms))
+      countsnap+=1
+      continue
+    if line.find('ITEM: ATOMS')==0:
+      read_atom_flag=True
+      continue
+    if read_atom_flag==True:
+      read_atoms+=1
+      parse=line.split()
+      if parse[0]!="":
+        coords[int(parse[0])-1][0]=int(parse[1])
+        coords[int(parse[0])-1][1]=float(parse[2])
+        coords[int(parse[0])-1][2]=float(parse[3])
+        coords[int(parse[0])-1][3]=float(parse[4])
+      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]))
+
+  print "\nDone converting to xyz!!\n"
+  fin.close()
+  fout.close()
+  return
+
+if __name__ == '__main__':
+
+    # if no input, print help and exit
+    if len(sys.argv) < 2:
+        print Info()
+        sys.exit(1)
+
+    inputfile=sys.argv[1]
+    outfile=sys.argv[2]
+
+    lmp2xyz(inputfile,outfile.split()[0])
+

tools/eff/radii.vmd

+# module radii.vmd
+# December, 2009 -(c)- Andres Jaramillo-Botero
+# Script to load variable changing radii onto a pEFF lammpstrj file
+
+# radii.vmd --
+#    Script to read and change electron radii in vmd dynamics traj
+#
+ 
+# openFile --
+#    Open the file and start looking for data
+#
+# Arguments:
+#    filename    Name of the file to read
+#
+# Result:
+#    infile      Handle to the opened file
+#
+# Side effects:
+#    A file in question is opened. Lines containing a * or# as
+#    the first non-blank character are considered comments. The
+#    first line without this is considered to be a line with the
+#    names of the columns.
+#
+proc openFile {filename} {
+    set infile [open $filename "r"]
+    return $infile
+}
+
+# readNames --
+#    Read the names of the columns
+#
+# Arguments:
+#    infile      Handle to the file
+#
+# Result:
+#    names       List of names, the number indicates the number of
+#                the snapshot frame
+#
+proc readNames {infile} {
+  #
+  # Skip the header - if any
+  #
+    set pos    0
+    while { [gets $infile line] >= 0 } {
+       if { [regexp {[ \t]*[*#]} $line] } {
+          incr pos
+       } else {
+          break
+       }
+    }
+    seek $infile 0 start
+    while { $pos > 0 } {
+       gets $infile line
+       incr pos -1
+    }
+
+  #
+  # Read the line with the column names
+  #
+    gets $infile line
+
+  # Force the line to be interpreted as a list
+    set nocols [llength $line]
+
+    return $line
+}
+
+# readData --
+#    Read the data per line
+#
+# Arguments:
+#    infile      Handle to the file
+#
+# Result:
+#    values      List of values, representing each column.
+#                A list of length zero indicates the end of file
+#
+proc readData {infile} {
+    while { [gets $infile values] == 0 } { ;# Just go on - skip empty lines }
+
+    set nocols [llength $values]
+
+    return $values
+}
+
+# readFile --
+#    Read the file and store the data in a (global) array
+#
+# Arguments:
+#    filename    Name of the file
+#
+# Result:
+#    None
+#
+# Side effects:
+#    Filled array, ready for display
+#
+proc readFile {filename} {
+    global data_array
+
+    set infile [openFile $filename]
+
+    set data_array(names) [readNames $infile]
+    set i 0
+    foreach name $data_array(names) {
+       set data_array($i) {}
+       incr i
+    }
+
+    while 1 {
+       set values [readData $infile]
+       if { [llength $values] > 0 } {
+          set i 0
+          foreach value $values {
+             lappend data_array($i) $value
+             incr i
+          }
+       } else {
+          break
+       }
+    }
+}
+
+# makeXYData --
+#    Make a list useable by frame-electron
+#
+# Arguments:
+#    xindex     Index of frame data
+#    yindex     Index of electron data
+#
+# Result:
+#    None
+#
+# Side effects:
+#    A dataset for changing the electron radii, per trajectory frame
+#
+proc makeXYData {xindex yindex} {
+    global data_array
+
+    set xydata {}
+    foreach x $data_array($xindex) y $data_array($yindex) {
+        lappend xydata $x $y
+    }
+    return $xydata
+}
+
+proc returnXYPair {x y} {
+   global data_array
+
+   return [list $data_array($x) $data_array($y)]
+}
+
+# do_radii --
+#    Changes the radii of electrons per trajectory frame
+#
+# Arguments:
+#    xindex     Index of frame data
+#    yindex     Index of electron data
+#
+# Result:
+#    prints the frame:atomID:radius
+#
+proc do_radii {args} {
+  global molid data_array
+
+  #set n [molinfo $molid get numatoms]
+  set f [molinfo $molid get frame]
+ 
+  set fr [expr {$f+1}]
+
+  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
+    $s set radius $nr
+    $s delete
+    #puts stderr "$fr $elec $nr"
+  }
+}
+
+# main --
+#    Main control flow
+#
+
+# Check input arguments
+#for {set i 0} {$i<[llength $argv]} {incr i}
+#  puts " - $i: [lindex $argv $i]"
+
+global data_array molid
+
+# Set input files manually
+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]
+
+# load electron radii information for trajectory
+#set datafile [lindex $::argv 1]
+
+set molid [mol new $xyz waitfor all]
+mol modstyle 0 [molinfo top] VDW 1.0 32.0
+
+puts "Starting ..."
+readFile $data
+puts "Read datafile and created array of radii ..."
+
+puts "Visualize trajectory" 
+trace variable vmd_frame($molid) w do_radii
+animate goto start
+do_radii
+  

tools/eff/setup.py

+from distutils.core import setup
+from distutils.extension import Extension
+from Cython.Distutils import build_ext
+
+# To build do: python setup.py build_ext --inplace
+
+ext_modules = [Extension("lmp2radii",["lmp2radii.pyx"])]
+setup(
+    name = 'lmp2radii',
+    cmdclass = {'build_ext': build_ext},
+    ext_modules = ext_modules
+)

tools/eff/structure_generators/Be-solid.pl

+#!/usr/bin/perl
+
+# Usage: Be-solid.pl <nx> <ny> <nz> > data.Be
+#
+# Generates a beryllium solid LAMMPS data file.
+
+$nx = shift(@ARGV);
+$ny = shift(@ARGV);
+$nz = shift(@ARGV);
+
+$La = 4.592;  # eFF optimized (4.33 expt)
+$Lc = 7.025;  # eFF optimized (6.78 expt)
+
+# This part changes for different lattices
+@xunit = (0, 0.5, 0, 0.5); 
+@yunit = (0, 0.5, 0.5 * (4/3), (1/3) * 0.5);
+@zunit = (0, 0, 0.5, 0.5);
+
+$Lx = $La;
+$Ly = $La * sqrt(3);
+$Lz = $Lc;
+
+$idx = 0;
+for ($x = 0; $x < $nx; $x++)
+{
+  for ($y = 0; $y < $ny; $y++)
+  {
+    for ($z = 0; $z < $nz; $z++)
+    {
+      for ($i = 0; $i <= $#xunit; $i++)
+      {
+        $xnuc[$idx] = ($x + .25) * $Lx + $xunit[$i] * $Lx;
+        $ynuc[$idx] = ($y + .25) * $Ly + $yunit[$i] * $Ly;
+        $znuc[$idx] = ($z + .25) * $Lz + $zunit[$i] * $Lz;
+        $idx++;
+      } 
+    }
+  }
+} 
+
+$numnuc = $idx;
+
+# Print length of supercell
+
+printf("Created by Be-solid.pl\n\n");
+printf("%d atoms\n",$numnuc+$numnuc*2+$numnuc*2);
+printf("2 atom types\n\n");
+printf("%f %f xlo xhi\n", 0, $Lx * $nx);
+printf("%f %f ylo yhi\n", 0, $Ly * $ny);
+printf("%f %f zlo zhi\n\n", 0, $Lz * $nz);
+printf("Masses\n\n");
+printf("1 9.012182\n");
+printf("2 1.000000\n\n");
+printf("Atoms\n\n");
+
+$j=0;
+# Print out the nuclei
+for ($i = 0; $i < $numnuc; $i++)
+{
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 1, 4.0, 0, 0.0, $xnuc[$i], $ynuc[$i], $znuc[$i]);
+}
+
+# Print out the core electrons
+$r_core = 0.5;
+for ($i = 0; $i < $numnuc; $i++)
+{
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2,  0.0, 1, $r_core,$xnuc[$i], $ynuc[$i], $znuc[$i]);
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 0.0, -1, $r_core,$xnuc[$i], $ynuc[$i], $znuc[$i]);
+}
+
+# Print out the valence electrons
+$r_valence = 2.5;
+for ($i = 0; $i < $numnuc; $i += 4)
+{
+  $x = &BoundX($xnuc[$i] + $Lx / 2.0);
+  $y = &BoundY($ynuc[$i] + $Lx / (2.0 * sqrt(3)));
+  $z = $znuc[$i];
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 0.0, 1, $r_valence,$x, $y, $z);
+  $x = &BoundX($xnuc[$i] + $Lx / 2.0);
+  $y = &BoundY($ynuc[$i] + $Lx / (2.0 * sqrt(3)));
+  $z = $znuc[$i];
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 0.0, -1, $r_valence,$x, $y, $z);
+
+  $x = &BoundX($xnuc[$i+1] + $Lx / 2.0);
+  $y = &BoundY($ynuc[$i+1] + $Lx / (2.0 * sqrt(3)));
+  $z = $znuc[$i+1];
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 0.0, 1, $r_valence, $x, $y, $z);
+  $x = &BoundX($xnuc[$i+1] + $Lx / 2.0);
+  $y = &BoundY($ynuc[$i+1] + $Lx / (2.0 * sqrt(3)));
+  $z = $znuc[$i+1];
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 0.0, -1, $r_valence,$x, $y, $z);
+
+  $x = &BoundX($xnuc[$i+2] - $Lx / 2.0);
+  $y = &BoundY($ynuc[$i+2] - $Lx / (2.0 * sqrt(3)));
+  $z = $znuc[$i+2];
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 0.0, 1, $r_valence,$x, $y, $z);
+  $x = &BoundX($xnuc[$i+2] - $Lx / 2.0);
+  $y = &BoundY($ynuc[$i+2] - $Lx / (2.0 * sqrt(3)));
+  $z = $znuc[$i+2];
+  printf("%i %i %f %f %f %f %f %f\n", $j+=1, 2, 0.0, -1, $r_valence,$x, $y, $z);
+
+  $x = &BoundX($xnuc[$i+3] - $Lx / 2.0);
+  $y = &BoundY($ynuc[$i+3] - $Lx / (2.0 * sqrt(3)));
+  $z = $znuc[$i+3];
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 0.0, 1, $r_valence, $x, $y, $z);
+  $x = &BoundX($xnuc[$i+3] - $Lx / 2.0);
+  $y = &BoundY($ynuc[$i+3] - $Lx / (2.0 * sqrt(3)));
+  $z = $znuc[$i+3];
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 0.0, -1, $r_valence,$x, $y, $z);
+}
+
+sub BoundX {
+  $x = $_[0];
+  if ($x > $Lx * $nx) {$x -= $Lx * $nx;}
+  if ($x < 0) {$x += $Lx * $nx;}
+  return $x;
+}
+
+sub BoundY {
+  $y = $_[0];
+  if ($y > $Ly * $ny) {$y -= $Ly * $ny;}
+  if ($y < 0) {$y += $Ly * $ny;}
+  return $y;
+}
+  
+sub BoundZ {
+  $z = $_[0];
+  if ($z > $Lz * $nz) {$z -= $Lz * $nz;}
+  if ($z < 0) {$z += $Lz * $nz;}
+  return $z;
+}

tools/eff/structure_generators/Diamond.pl

+#!/usr/bin/perl
+
+# Usage: diamond <nx> <ny> <nz> > data.diamond 
+#
+# Generates diamond A4 structure.
+
+$nx = shift(@ARGV);
+$ny = shift(@ARGV);
+$nz = shift(@ARGV);
+#$scale = shift(@ARGV);
+
+$scale = 1.045;
+$L = 7.33461;       # eFF optimized (6.740 expt)
+$re_core  = 0.328;        # core electron
+$re_sigma = 1.559;        # sigma electrons
+
+#$r_threshold = $scale * 1.6 / 0.529;          # threshold for sigma bond length
+$r_threshold = $scale * 1.6 / 0.5 ;          # threshold for sigma bond length
+
+@xunit = (0, 0, 0.5, 0.5, 0.25, 0.75, 0.25, 0.75);
+@yunit = (0, 0.5, 0, 0.5, 0.25, 0.75, 0.75, 0.25);
+@zunit = (0, 0.5, 0.5, 0, 0.25, 0.25, 0.75, 0.75);
+
+$idx = 0;
+for ($x = 0; $x < $nx; $x++)
+{
+  for ($y = 0; $y < $ny; $y++)
+  {
+    for ($z = 0; $z < $nz; $z++)
+    {
+      for ($i = 0; $i <= $#xunit; $i++)
+      {
+        $xnuc[$idx] = $x * $L + $xunit[$i] * $L;
+        $ynuc[$idx] = $y * $L + $yunit[$i] * $L;
+        $znuc[$idx] = $z * $L + $zunit[$i] * $L;
+        $idx++;
+      } 
+    }
+  }
+} 
+
+$numnuc = $idx;
+$xx = $L * $nx;
+$yy = $L * $ny;
+$zz = $L * $nz;
+
+# Print length of supercell
+
+printf("Created with Diamond.pl\n\n");
+printf("%d atoms\n",$numnuc*7);
+printf("2 atom types\n\n");
+printf("%f %f xlo xhi\n", 0, $L * $nx);
+printf("%f %f ylo yhi\n", 0, $L * $ny);
+printf("%f %f zlo zhi\n\n", 0, $L * $nz);
+printf("Masses\n\n");
+printf("1 12.01070\n");
+printf("2 1.000000\n\n");
+printf("Atoms\n\n");
+
+# Print out the nuclei and the core electrons
+for ($i = 0; $i < $numnuc; $i++)
+{
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 1, 6.0, 0.0, 0.0, $xnuc[$i], $ynuc[$i], $znuc[$i]);
+}
+
+for ($i = 0; $i < $numnuc; $i++)
+{
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 0.0, 1, $re_core,$xnuc[$i], $ynuc[$i], $znuc[$i]);
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 0.0,-1, $re_core,$xnuc[$i], $ynuc[$i], $znuc[$i]);
+}
+
+# Print out sigma electrons
+
+$LLx = $L * $nx;
+$LLy = $L * $ny;
+$LLz = $L * $nz;
+
+$k=$j;
+for ($i = 0; $i < $numnuc; $i++)
+{
+  for ($j = 0; $j < $i; $j++)
+  {
+    $x = $xnuc[$j] - $xnuc[$i];
+    $y = $ynuc[$j] - $ynuc[$i];
+    $z = $znuc[$j] - $znuc[$i];
+
+    # Minimum image convention
+
+    if ($x > $LLx/2) {$x -= $LLx;}
+    if ($y > $LLy/2) {$y -= $LLy;}
+    if ($z > $LLz/2) {$z -= $LLz;}
+    
+    if ($x < -$LLx/2) {$x += $LLx;}
+    if ($y < -$LLy/2) {$y += $LLy;}
+    if ($z < -$LLz/2) {$z += $LLz;}
+
+    $r = sqrt($x * $x + $y * $y + $z * $z);
+    if ($r   < $r_threshold)
+    {
+      $bond_x = $xnuc[$i] + $x / 2;
+      $bond_y = $ynuc[$i] + $y / 2;
+      $bond_z = $znuc[$i] + $z / 2;
+
+      # Minimum image convention
+
+      if ($bond_x > $LLx) {$bond_x -= $LLx};
+      if ($bond_y > $LLy) {$bond_y -= $LLy};
+      if ($bond_z > $LLz) {$bond_z -= $LLz};
+
+      if ($bond_x < 0)   {$bond_x += $LLx};
+      if ($bond_y < 0)   {$bond_y += $LLy};
+      if ($bond_z < 0)   {$bond_z += $LLz};
+
+      printf("%i %i %f %i %f %f %f %f\n", $k+=1, 2, 0.0, 1, $re_sigma,$bond_x, $bond_y, $bond_z);
+      printf("%i %i %f %i %f %f %f %f\n", $k+=1, 2, 0.0,-1, $re_sigma,$bond_x, $bond_y, $bond_z);
+    }
+  }
+}
+

tools/eff/structure_generators/Li-hydride.pl

+#!/usr/bin/perl
+
+# Usage: ./Li-hydride.pl <nx> <ny> <nz> > data.Li-hydride
+#
+# Generates a lithium hydride solid.
+
+$nx = shift(@ARGV);
+$ny = shift(@ARGV);
+$nz = shift(@ARGV);
+
+$L = 6.785;  # eFF optimized (7.72 expt)
+
+# This part changes for different lattices
+@xunit = (0, 0.5, 0, 0.5, 0.5, 0, 0, 0.5);
+@yunit = (0, 0.5, 0.5, 0, 0, 0.5, 0, 0.5);
+@zunit = (0, 0, 0.5, 0.5, 0, 0, 0.5, 0.5);
+
+$r_elec = 0.7;
+$r2_elec = 2.2;
+
+$Lx = $L;
+$Ly = $L;
+$Lz = $L;
+
+$idx = 0;
+for ($x = 0; $x < $nx; $x++)
+{
+  for ($y = 0; $y < $ny; $y++)
+  {
+    for ($z = 0; $z < $nz; $z++)
+    {
+      for ($i = 0; $i <= $#xunit; $i++)
+      {
+        $xnuc[$idx] = $x * $Lx + $xunit[$i] * $L + 0.5 * $L;
+        $ynuc[$idx] = $y * $Ly + $yunit[$i] * $L + 0.5 * $L;
+        $znuc[$idx] = $z * $Lz + $zunit[$i] * $L + 0.5 * $L;
+        $idx++;
+      } 
+    }
+  }
+} 
+
+$numnuc = $idx;
+
+# Print length of supercell
+
+printf("Created with Li-hydride.pl\n\n");
+printf("%d atoms\n",$numnuc+$numnuc+$numnuc);
+printf("3 atom types\n\n");
+printf("%f %f xlo xhi\n", 0, $Lx * $nx);
+printf("%f %f ylo yhi\n", 0, $Ly * $ny);
+printf("%f %f zlo zhi\n\n", 0, $Lz * $nz);
+printf("Masses\n\n");
+printf("1 6.941000\n");
+printf("2 1.007940\n");
+printf("3 1.000000\n\n");
+printf("Atoms\n\n");
+
+$j = 0;
+# Print out the nuclei
+for ($i = 0; $i < $numnuc; $i += 8)
+{
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 1, 3.0, 0, 0.0,$xnuc[$i], $ynuc[$i], $znuc[$i]);
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 1, 3.0, 0, 0.0,$xnuc[$i+1], $ynuc[$i+1], $znuc[$i+1]);
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 1, 3.0, 0, 0.0,$xnuc[$i+2], $ynuc[$i+2], $znuc[$i+2]);
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 1, 3.0, 0, 0.0,$xnuc[$i+3], $ynuc[$i+3], $znuc[$i+3]);
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 1.0, 0, 0.0,$xnuc[$i+4], $ynuc[$i+4], $znuc[$i+4]);
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 1.0, 0, 0.0,$xnuc[$i+5], $ynuc[$i+5], $znuc[$i+5]);
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 1.0, 0, 0.0,$xnuc[$i+6], $ynuc[$i+6], $znuc[$i+6]);
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 1.0, 0, 0.0,$xnuc[$i+7], $ynuc[$i+7], $znuc[$i+7]);
+}
+
+# Print out the electrons
+for ($i = 0; $i < $numnuc; $i += 8)
+{
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 3, 0.0, 1, $r_elec,$xnuc[$i  ], $ynuc[$i  ], $znuc[$i  ]);
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 3, 0.0, -1, $r_elec,$xnuc[$i  ], $ynuc[$i  ], $znuc[$i  ]);
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 3, 0.0, 1, $r_elec,$xnuc[$i+1], $ynuc[$i+1], $znuc[$i+1]);
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 3, 0.0, -1, $r_elec,$xnuc[$i+1], $ynuc[$i+1], $znuc[$i+1]);
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 3, 0.0, 1, $r_elec,$xnuc[$i+2], $ynuc[$i+2], $znuc[$i+2]);
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 3, 0.0, -1, $r_elec,$xnuc[$i+2], $ynuc[$i+2], $znuc[$i+2]);
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 3, 0.0, 1, $r_elec,$xnuc[$i+3], $ynuc[$i+3], $znuc[$i+3]);
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 3, 0.0, -1, $r_elec,$xnuc[$i+3], $ynuc[$i+3], $znuc[$i+3]);
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 3, 0.0, 1, $r_elec,$xnuc[$i+4], $ynuc[$i+4], $znuc[$i+4]);
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 3, 0.0, -1, $r_elec,$xnuc[$i+4], $ynuc[$i+4], $znuc[$i+4]);
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 3, 0.0, 1, $r_elec,$xnuc[$i+5], $ynuc[$i+5], $znuc[$i+5]);
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 3, 0.0, -1, $r_elec,$xnuc[$i+5], $ynuc[$i+5], $znuc[$i+5]);
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 3, 0.0, 1, $r_elec,$xnuc[$i+6], $ynuc[$i+6], $znuc[$i+6]);
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 3, 0.0, -1, $r_elec,$xnuc[$i+6], $ynuc[$i+6], $znuc[$i+6]);
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 3, 0.0, 1, $r_elec,$xnuc[$i+7], $ynuc[$i+7], $znuc[$i+7]);
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 3, 0.0, 1, $r_elec,$xnuc[$i+7], $ynuc[$i+7], $znuc[$i+7]);
+}
+

tools/eff/structure_generators/Li-solid.pl

+#!/usr/bin/perl
+
+# Usage: Li-solid.pl <nx> <ny> <nz> > data.Li
+#
+# Generates a lithium solid LAMMPS data input file
+# with FCC nuclear positions and interstitial electron positions.
+
+$nx = shift(@ARGV);
+$ny = shift(@ARGV);
+$nz = shift(@ARGV);
+
+$L = 8.352;  # eFF optimized (8.32 expt)
+
+# This part changes for different lattices
+@xunit = (0, 0.5, 0, 0.5, 0.5, 0, 0, 0.5);
+@yunit = (0, 0.5, 0.5, 0, 0, 0.5, 0, 0.5);
+@zunit = (0, 0, 0.5, 0.5, 0, 0, 0.5, 0.5);
+
+$r_elec = 0.7;
+$r_elec2 = 3.0;
+
+$Lx = $L;
+$Ly = $L;
+$Lz = $L;
+
+$idx = 0;
+for ($x = 0; $x < $nx; $x++)
+{
+ for ($y = 0; $y < $ny; $y++)
+ {
+   for ($z = 0; $z < $nz; $z++)
+   {
+     for ($i = 0; $i <= $#xunit; $i++)
+     {
+       $xnuc[$idx] = $x * $Lx + $xunit[$i] * $L + 0.5 * $L;
+       $ynuc[$idx] = $y * $Ly + $yunit[$i] * $L + 0.5 * $L;
+       $znuc[$idx] = $z * $Lz + $zunit[$i] * $L + 0.5 * $L;
+       $idx++;
+     } 
+   }
+ }
+} 
+
+$numnuc = $idx;
+
+# Print length of supercell
+
+printf("Created by AJB\n\n");
+printf("%d atoms\n",2*$numnuc);
+printf("2 atom types\n\n");
+printf("%f %f xlo xhi\n", 0, $Lx * $nx);
+printf("%f %f ylo yhi\n", 0, $Ly * $ny);
+printf("%f %f zlo zhi\n\n", 0, $Lz * $nz);
+printf("Masses\n\n");
+printf("1 6.941000\n");
+printf("2 1.000000\n\n");
+printf("Atoms\n\n");
+
+$j = 0;
+# Print out the nuclei
+for ($i = 0; $i < $numnuc; $i += 8)
+{
+ printf("%i %i %f %i %f %f %f %f\n", $j+=1,1,3.0, 0, 0.0,$xnuc[$i], $ynuc[$i], $znuc[$i]);
+ printf("%i %i %f %i %f %f %f %f\n", $j+=1,1,3.0, 0, 0.0,$xnuc[$i+1], $ynuc[$i+1], $znuc[$i+1]);
+ printf("%i %i %f %i %f %f %f %f\n", $j+=1,1,3.0, 0, 0.0,$xnuc[$i+2], $ynuc[$i+2], $znuc[$i+2]);
+ printf("%i %i %f %i %f %f %f %f\n", $j+=1,1,3.0, 0, 0.0,$xnuc[$i+3], $ynuc[$i+3], $znuc[$i+3]);
+}
+
+# Print out the core electrons
+for ($i = 0; $i < $numnuc; $i += 8)
+{
+ printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 0.0, 1, $r_elec,$xnuc[$i  ], $ynuc[$i  ], $znuc[$i  ]);
+ printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 0.0, -1, $r_elec,$xnuc[$i  ], $ynuc[$i  ], $znuc[$i  ]);
+ printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 0.0, 1, $r_elec,$xnuc[$i+1], $ynuc[$i+1], $znuc[$i+1]);
+ printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 0.0, -1, $r_elec,$xnuc[$i+1], $ynuc[$i+1], $znuc[$i+1]);
+ printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 0.0, 1, $r_elec,$xnuc[$i+2], $ynuc[$i+2], $znuc[$i+2]);
+ printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 0.0, -1, $r_elec,$xnuc[$i+2], $ynuc[$i+2], $znuc[$i+2]);
+ printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 0.0, 1, $r_elec,$xnuc[$i+3], $ynuc[$i+3], $znuc[$i+3]);
+ printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 0.0, -1, $r_elec,$xnuc[$i+3], $ynuc[$i+3], $znuc[$i+3]);
+}
+
+# Print out the valence electrons
+for ($i = 0; $i < $numnuc; $i += 8)
+{
+ if (rand() < .5) {$spin = 1;} else {$spin = -1;}  
+ printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 0.0, $spin, $r_elec2,$xnuc[$i+4], $ynuc[$i+4], $znuc[$i+4]);
+ printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 0.0, -$spin, $r_elec2,$xnuc[$i+5], $ynuc[$i+5], $znuc[$i+5]);
+ if (rand() < .5) {$spin = 1;} else {$spin = -1;}  
+ printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 0.0, $spin, $r_elec2,$xnuc[$i+6], $ynuc[$i+6], $znuc[$i+6]);
+ printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 0.0, -$spin, $r_elec2,$xnuc[$i+7], $ynuc[$i+7], $znuc[$i+7]);
+}
+

tools/eff/structure_generators/Uniform-electron-gas.pl

+#!/usr/bin/perl
+
+# Usage: ./Uniform-electron-gas.pl <nx> <ny> <nz> <rs> > data.uniform-e-gas
+#
+# Generates a uniform electron gas using electrons on an NaCl lattice.
+
+$nx = shift(@ARGV);
+$ny = shift(@ARGV);
+$nz = shift(@ARGV);
+$rs = shift(@ARGV);
+
+# This part changes for different lattices
+@xunit = (0, 0.5, 0.5, 0, 0, 0.5, 0.5, 0);
+@yunit = (0, 0, 0.5, 0.5, 0.5, 0.5, 0, 0);
+@zunit = (0, 0, 0, 0, 0.5, 0.5, 0.5, 0.5);
+$volume = 1;
+$elecs_per_cell = 8;
+
+$L = $rs * ((4/3) * 3.14159265 * $elecs_per_cell / $volume) ** (1/3);  # length of unit cell
+$r_elec = $rs * sqrt(1.5 / 1.1050);
+
+$Lx = $L;
+$Ly = $L;
+$Lz = $L;
+
+$idx = 0;
+for ($x = 0; $x < $nx; $x++)
+{
+  for ($y = 0; $y < $ny; $y++)
+  {
+    for ($z = 0; $z < $nz; $z++)
+    {
+      for ($i = 0; $i <= $#xunit; $i++)
+      {
+        $xnuc[$idx] = $x * $Lx + $xunit[$i] * $L;
+        $ynuc[$idx] = $y * $Ly + $yunit[$i] * $L;
+        $znuc[$idx] = $z * $Lz + $zunit[$i] * $L;
+        $idx++;
+      } 
+    }
+  }
+} 
+
+$numnuc = $idx;
+
+# Print length of supercell
+
+printf("Created with Uniform-electron-gas.pl\n\n");
+printf("%d atoms\n",$numnuc);
+printf("1 atom types\n\n");
+printf("%f %f xlo xhi\n", 0, $Lx * $nx);
+printf("%f %f ylo yhi\n", 0, $Ly * $ny);
+printf("%f %f zlo zhi\n\n", 0, $Lz * $nz);
+printf("Masses\n\n");
+printf("1 1.000000\n\n");
+printf("Atoms\n\n");
+
+$j = 0;
+
+# Print out the electrons
+for ($i = 0; $i < $numnuc; $i += 2)
+{
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 0.0, 1, 1, $r_elec, $xnuc[$i], $ynuc[$i], $znuc[$i]);
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 0.0, 1, -1, $r_elec, $xnuc[$i+1], $ynuc[$i+1], $znuc[$i+1]);
+}
+
+printf("\n");

tools/eff/structure_generators/h2.pl

+#!/usr/bin/perl
+
+# Usage: h2.pl <nx> <ny> <nz> <rs> > data.h2
+#
+# Generates rectangular lattice of hydrogen atoms, # in each direction = nx,
+# with Wigner radius rs.
+#    ...
+# Shifts H atoms in alternating directions so they form H2 molecule
+# starting structures.
+
+$nx = shift(@ARGV);
+$ny = shift(@ARGV);
+$nz = shift(@ARGV);
+$rs = shift(@ARGV);
+
+$L  = 1.6119919540164693 * $rs;  # length of unit cell ((4/3 pi)^(1/3)
+$re = 1.823572;                  # electron radius
+
+$dshift = 0.5;
+
+$idx = 0;
+for ($x = 0; $x < $nx; $x++)
+{
+  for ($y = 0; $y < $ny; $y++)
+  {
+    $dsign = 1; 
+    for ($z = 0; $z < $nz; $z++)
+    {
+      $xnuc[$idx] = $x * $L + 0.5 * $L; 
+      $ynuc[$idx] = $y * $L + 0.5 * $L;
+      $znuc[$idx] = $z * $L + 0.5 * $L + $dshift * $dsign;
+      $dsign = -$dsign;
+      $idx++;
+    }
+  }
+} 
+
+$numnuc = $idx;
+
+# Print length of supercell
+
+printf("Created with h2.pl\n\n");
+printf("%d atoms\n",$numnuc*2);
+printf("2 atom types\n\n");
+printf("%f %f xlo xhi\n", 0, $L * $nx);
+printf("%f %f ylo yhi\n", 0, $L * $ny);
+printf("%f %f zlo zhi\n\n", 0, $L * $nz);
+printf("Masses\n\n");
+printf("1 1.007940\n");
+printf("2 1.000000\n\n");
+printf("Atoms\n\n");
+
+$j=0;
+
+# Print out the nuclei and the core electrons
+for ($i = 0; $i < $numnuc; $i++)
+{
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 1, 1.0, 0, 0.0, $xnuc[$i], $ynuc[$i], $znuc[$i]);
+}
+
+$spin = 1;
+for ($i = 0; $i < $numnuc; $i++)
+{
+  if ($spin == 1) {$spin = -1;} else {$spin = 1;}
+  printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2,  0.0, $spin, $re, $xnuc[$i], $ynuc[$i], $znuc[$i]);
+}
+