Merge branch 'master' of github.com:lkwagner/autogenv2

Author: lkwagner

Crystal.py

@@ -18,7 +18,7 @@ def __init__(self):
     #Electron model
     self.spin_polarized=True    
     self.xml_name="BFD_Library.xml"
-    self.functional={'exchange':'PBE','correlation':'PBE','hybrid':0}
+    self.functional={'exchange':'PBE','correlation':'PBE','hybrid':0,'predefined':None}
     self.total_spin=0
 
     #Numerical convergence parameters
@@ -27,7 +27,7 @@ def __init__(self):
     self.kmesh=[8,8,8]
     self.gmesh=16
     self.tolinteg=[8,8,8,8,18]
-    self.dftgrid='XLGRID'
+    self.dftgrid=''
 
     #Memory
     self.biposize=100000000
@@ -41,7 +41,9 @@ def __init__(self):
     self.levshift=[]
     self.broyden=[0.01,60,15]
     self.diis=False
-    self.smear=0.0001
+    self.smear=0.0
+    self.supercell= [[1.,0.,0.],[0.,1.,0.],[0.,0.,1.]]  
+    
 
     # Use the new crystal2qmc script. This should change soon!
     self.cryapi=True
@@ -55,7 +57,6 @@ def set_struct_fromcif(self,cifstr,primitive=True):
     self.primitive=primitive
     self.cif=cifstr
     self.struct=CifParser.from_string(self.cif).get_structures(primitive=self.primitive)[0].as_dict()
-    self.supercell= [[1.,0.,0.],[0.,1.,0.],[0.,0.,1.]]  
     self.boundary="3d"
   #-----------------------------------------------
 
@@ -75,7 +76,7 @@ def set_options(self, d):
       selfdict[k]=d[k]
 
   #-----------------------------------------------
-  def crystal_input(self):
+  def crystal_input(self,section4=[]):
 
     assert self.struct is not None,'Need to set "struct" first.'
     geomlines=self.geom()
@@ -99,15 +100,19 @@ def crystal_input(self):
     outlines+=["DFT"]
     if self.spin_polarized:
       outlines+=["SPIN"]
-    outlines += [ 
-      "EXCHANGE",
-      self.functional['exchange'],
-      "CORRELAT",
-      self.functional['correlation'],
-      "HYBRID", 
-      str(self.functional['hybrid']),
-      self.dftgrid,
-      "END",
+    if self.functional['predefined']!=None:
+      outlines+=[self.functional['predefined']]
+    else:
+      outlines += [ 
+        "EXCHANGE",
+        self.functional['exchange'],
+        "CORRELAT",
+        self.functional['correlation'],
+        "HYBRID", 
+        str(self.functional['hybrid'])]
+    if self.dftgrid!="":
+      outlines+=[self.dftgrid]
+    outlines+=["END",
       "SCFDIR",
       "BIPOSIZE",
       str(self.biposize),
@@ -121,10 +126,10 @@ def crystal_input(self):
       ' '.join(map(str,self.tolinteg)),
       "MAXCYCLE",
       str(self.maxcycle),
-      "SMEAR",
-      str(self.smear),
       "SAVEWF"
     ]
+    if self.smear > 0:
+      outlines+=["SMEAR",str(self.smear)]
     if self.spin_polarized:
       outlines+=['SPINLOCK',str(self.total_spin)+" 200"]
 
@@ -134,7 +139,7 @@ def crystal_input(self):
       outlines+=["LEVSHIFT"," ".join(map(str,self.levshift))]
     else:
       outlines+=["BROYDEN"," ".join(map(str,self.broyden))]
-
+    outlines+=section4
     if self.restart:
       outlines+=["GUESSP"]
     outlines+=["END"]
@@ -262,19 +267,13 @@ def basis_section(self):
       elements.add(nm)
     basislines=[]
     elements = sorted(list(elements)) # Standardize ordering.
-    if self.boundary=='3d':
-      for e in elements:
-        basislines+=self.generate_basis_3d(e)
-    elif self.boundary=='0d':
-      for e in elements:
-        basislines+=self.generate_basis_0d(e)
-    else:
-      raise AssertionError("Invalid boundary value")
+    for e in elements:
+      basislines+=self.generate_basis(e)
     return basislines
 
 
 ########################################################
-  def generate_basis_3d(self,symbol):
+  def generate_basis(self,symbol):
     """
     Author: "Kittithat (Mick) Krongchon" <[email protected]> and Lucas K. Wagner
     Returns a string containing the basis section.  It is modified according to a simple recipe:
@@ -360,75 +359,6 @@ def generate_basis_3d(self,symbol):
             ret+=["0 %i %i %g 1"%(basis_index[angular],1,0.0),line]
             ncontract+=1
 
-    return ["%i %i"%(Element(symbol).number+200,ncontract)] +\
-            self.pseudopotential_section(symbol) +\
-            ret
-
-########################################################
-  def generate_basis_0d(self,symbol):
-    """
-    Author: "Kittithat (Mick) Krongchon" <[email protected]> and Lucas K. Wagner
-    Returns a string containing the basis section.  It is modified according to a simple recipe:
-    Args:
-        symbol (str): The symbol of the element to be specified in the
-            D12 file.
-        xml_name (str): The name of the XML pseudopotential and basis
-            set database.
-    Returns:
-        str: The pseudopotential and basis section.
-    """
-    
-    maxorb=3
-    basis_name="vtz"
-    maxcharge={"s":2,"p":6,"d":10,"f":15} #max charge on a basis function
-    max_tot_charge={"s":2,"p":6,"d":10,"f":15} #max total charge in an angular momentum
-    
-    basis_index={"s":0,"p":2,"d":3,"f":4}
-    transition_metals=["Sc","Ti","V","Cr","Mn","Fe","Co","Ni","Cu","Zn"]
-    if symbol in transition_metals:
-      maxorb=4
-      max_tot_charge['s']=4
-    charge_sum={"s":0,"p":0,"d":0,"f":0}
-    tree = ElementTree()
-    tree.parse(self.xml_name)
-    element = tree.find('./Pseudopotential[@symbol="{}"]'.format(symbol))
-    atom_charge = int(element.find('./Effective_core_charge').text)
-    if symbol in self.initial_charges.keys():
-      atom_charge-=self.initial_charges[symbol]
-    basis_path = './Basis-set[@name="{}"]/Contraction'.format(basis_name)
-    found_orbitals = []
-    totcharge=0
-    ret=[]
-    ncontract=0
-    for contraction in element.findall(basis_path):
-        angular = contraction.get('Angular_momentum')
-
-        #Figure out which coefficients to print out based on the minimal exponent
-        nterms = 0
-        basis_part=[]
-        for basis_term in contraction.findall('./Basis-term'):
-            exp = basis_term.get('Exp')
-            coeff = basis_term.get('Coeff')
-            basis_part += ['  {} {}'.format(exp, coeff)]
-            nterms+=1
-        #now write the header 
-        if nterms > 0 and angular in basis_index.keys():
-          found_orbitals.append(angular)          
-          charge=min(atom_charge-totcharge,maxcharge[angular])
-          if charge_sum[angular] >= max_tot_charge[angular]:
-            charge=0
-          #put in a special case for transition metals:
-          #depopulate the 4s if the atom is charged
-          if symbol in transition_metals and symbol in self.initial_charges.keys() \
-                  and self.initial_charges[symbol] > 0 and found_orbitals.count(angular) > 1 \
-                  and angular=="s":
-            charge=0
-          totcharge+=charge
-          charge_sum[angular]+=charge
-          ret+=["0 %i %i %g 1"%(basis_index[angular],nterms,charge)] + basis_part
-          ncontract+=1
-
-
     return ["%i %i"%(Element(symbol).number+200,ncontract)] +\
             self.pseudopotential_section(symbol) +\
             ret
@@ -495,7 +425,6 @@ def collect(self,outfilename):
       lines = f.readlines()
       for li,line in enumerate(lines):
         if 'SCF ENDED' in line:
-          print(line)
           self.out['total_energy']=float(line.split()[8])    
 
         elif 'TOTAL ATOMIC SPINS' in line:
@@ -505,7 +434,8 @@ def collect(self,outfilename):
             moms += map(float,lines[li+shift].split())
             shift += 1
           self.out['mag_moments']=moms
-      print(self.out)
+      self.out['etots'] = [float(line.split()[3]) for line in lines if "DETOT" in line]
+      
       self.completed=True
     else:
       # Just to be sure/clear...
@@ -523,36 +453,36 @@ def write_summary(self):
   def check_outputfile(self,outfilename,acceptable_scf=10.0):
     """ Check output file. 
 
-    Current return values:
+    Return values:
     no_record, not_started, ok, too_many_cycles, finished (fall-back),
     scf_fail, not_enough_decrease, divergence, not_finished
     """
     if os.path.isfile(outfilename):
-      outf = open(outfilename,'r')
+      outf = open(outfilename,'r',errors='ignore')
     else:
       return "not_started"
 
-    outlines = outf.read().split('\n')
+    outlines = outf.readlines()
     reslines = [line for line in outlines if "ENDED" in line]
 
     if len(reslines) > 0:
       if "CONVERGENCE" in reslines[0]:
         return "ok"
       elif "TOO MANY CYCLES" in reslines[0]:
-        print("CrystalRunner: Too many cycles.")
+        #print("CrystalRunner: Too many cycles.")
         return "too_many_cycles"
       else: # What else can happen?
-        print("CrystalReader: Finished, but unknown state.")
+        #print("CrystalReader: Finished, but unknown state.")
         return "finished"
 
     detots = [float(line.split()[5]) for line in outlines if "DETOT" in line]
     if len(detots) == 0:
-      print("CrystalRunner: Last run completed no cycles.")
+      #print("CrystalRunner: Last run completed no cycles.")
       return "scf_fail"
 
     detots_net = sum(detots[1:])
     if detots_net > acceptable_scf:
-      print("CrystalRunner: Last run performed poorly.")
+      #print("CrystalRunner: Last run performed poorly.")
       return "not_enough_decrease"
 
     etots = [float(line.split()[3]) for line in outlines if "DETOT" in line]

crystal2qmc.py

@@ -373,6 +373,47 @@ def write_slater(basis,eigsys,kpt,base="qwalk",kfmt='coord'):
     outf.write("\n".join(outlines))
   return outlines # Might be confusing.
 
+###############################################################################
+def write_orbplot(basis,eigsys,kpt,base="qwalk",kfmt='coord'):
+  if kfmt == 'int': kbase = base + '_' + "{}".format(eigsys['kpt_index'][kpt])
+  else:             kbase = base + '_' + "{}{}{}".format(*kpt)
+  ntot = basis['ntot']
+  nmo  = basis['nmo']
+  nup  = eigsys['nup']
+  ndn  = eigsys['ndn']
+  uporbs = np.arange(nup)+1
+  dnorbs = np.arange(ndn)+1
+  if eigsys['nspin'] > 1:
+    dnorbs += nmo
+  if eigsys['ikpt_iscmpx'][kpt]: orbstr = "corbitals"
+  else:                          orbstr = "orbitals"
+  uporblines = ["{:5d}".format(orb) for orb in uporbs]
+  width = 10
+  for i in reversed(range(width,len(uporblines),width)):
+    uporblines.insert(i,"\n ")
+  dnorblines = ["{:5d}".format(orb) for orb in dnorbs]
+  for i in reversed(range(width,len(dnorblines),width)):
+    dnorblines.insert(i,"\n ")
+  outlines_prefix = [
+      "method { ",
+      "plot",
+      "{0} {{".format(orbstr),
+      "cutoff_mo",
+      "  magnify 1",
+      "  nmo {0}".format(dnorbs[-1]),
+      "  orbfile {0}.orb".format(kbase),
+      "  include {0}.basis".format(base),
+      "  centers { useglobal }",
+      "}",
+      "plotorbitals {",
+    ]
+  outlines_postfix= ["}","}","include "+kbase+".sys"]
+  with open(kbase+".up.plot",'w') as outf:
+    outf.write("\n".join(outlines_prefix+ [" ".join(uporblines)] + outlines_postfix))
+  with open(kbase+".dn.plot",'w') as outf:
+    outf.write("\n".join(outlines_prefix+ [" ".join(dnorblines)] + outlines_postfix))
+  
+
 ###############################################################################
 # f orbital normalizations are from 
 # <http://winter.group.shef.ac.uk/orbitron/AOs/4f/equations.html>
@@ -727,6 +768,7 @@ def convert_crystal(
   for kpt in eigsys['kpt_coords']:
     if eigsys['ikpt_iscmpx'][kpt] and kset=='real': continue
     write_slater(basis,eigsys,kpt,base,kfmt)
+    write_orbplot(basis,eigsys,kpt,base,kfmt)
     normalize_eigvec(eigsys,basis,kpt)
     write_orb(eigsys,basis,ions,kpt,base,kfmt)
     write_sys(lat_parm,basis,eigsys,pseudo,ions,kpt,base,kfmt)

example/make_crystal_plan.py

@@ -9,7 +9,6 @@
 
 dft_opts={
     'xml_name':os.getcwd()+'/../BFD_Library.xml',
-    'basis_params':[0.2,0,3],
     'cutoff':0.0,
     'dftgrid':'LGRID',
     'spin_polarized':False
@@ -29,7 +28,7 @@
 
 pbc_opts={
     'xml_name':os.getcwd()+'/../BFD_Library.xml',
-    'basis_params':[0.4,1,3],
+    'basis_params':[0.2,2,3],
     'cutoff':0.2,
     'dftgrid':'LGRID',
     'spin_polarized':False,