gamess analyzer, trimer structure

git-svn-id: https://subversion.assembla.com/svn/qmcdev/trunk@6562 e5b18d87-469d-4833-9cc0-8cdfa06e9491
This commit is contained in:
Jaron Krogel 2015-08-24 13:54:43 +00:00
parent dea5d759ca
commit 6fed27af80
6 changed files with 409 additions and 41 deletions

View File

@ -105,7 +105,7 @@ class TextFile(DevBase):
if pos!=-1:
return self.mm.seek(pos,0)
else:
return None
return -1
#end if
else:
return self.mm.seek(pos,whence)

View File

@ -33,6 +33,26 @@ def assign_value(host,dest,file,string):
class GamessAnalyzer(SimulationAnalyzer):
lset_full = 'spdfg'
lxyz = obj(
s = set('s'.split()),
p = set('x y z'.split()),
d = set('xx yy zz xy xz yz'.split()),
f = set('xxx yyy zzz xxy xxz yyx yyz zzx zzy xyz'.split()),
g = set('xxxx yyyy zzzz xxxy xxxz yyyx yyyz zzzx zzzy xxyy xxzz yyzz xxyz yyxz zzxy'.split()),
)
lxyz_reverse = obj()
for l,xyzs in lxyz.iteritems():
for xyz in xyzs:
lxyz_reverse[xyz] = l
#end for
#end for
def __init__(self,arg0=None,prefix=None,analyze=False,exit=False,**outfilenames):
self.info = obj(
exit = exit,
@ -111,27 +131,11 @@ class GamessAnalyzer(SimulationAnalyzer):
def analyze_log(self):
# read the log file
log = self.get_output('output')
# try to get the energy components
energy = obj()
try:
if log!=None and log.seek('ENERGY COMPONENTS',0)!=-1:
for n in xrange(18):
line = log.readline()
if '=' in line and 'ENERGY' in line:
nameline,value = line.split('=')
tokens = nameline.lower().split()
name = ''
for token in tokens:
if token!='energy':
name += token.replace('-','_')+'_'
#end if
#end for
name = name[:-1]
value = float(value.strip())
energy[name]=value
#end if
#end for
#end if
self.read_energy_components(log,energy)
except:
if self.info.exit:
self.error('log file analysis failed (energy components)')
@ -140,23 +144,279 @@ class GamessAnalyzer(SimulationAnalyzer):
if len(energy)>0:
self.energy = energy
#end if
# try to get the orbital count from the log file
counts = obj()
try:
if log.seek('TOTAL NUMBER OF MOS IN VARIATION SPACE',0)!=-1:
if log.seek('TOTAL NUMBER OF BASIS SET SHELLS',0)!=-1:
counts.shells = int(log.readtokens()[-1])
#end if
if log.seek('NUMBER OF CARTESIAN ATOMIC ORBITALS',0)!=-1:
counts.cao = int(log.readtokens()[-1])
#end if
if log.seek('TOTAL NUMBER OF MOS IN VARIATION SPACE',1)!=-1:
counts.mos = int(log.readtokens()[-1])
#end if
except:
if self.info.exit:
self.error('log file analysis failed (mo count)')
self.error('log file analysis failed (counts)')
#end if
#end try
if len(counts)>0:
self.counts = counts
#end if
# try to get the up/down orbitals
if 'counts' in self and False: # don't read orbitals, large
orbitals = obj()
try:
self.read_orbitals(log,orbitals,'up' ,'-- ALPHA SET --')
self.read_orbitals(log,orbitals,'down','-- BETA SET --' )
except:
if self.info.exit:
self.error('log file analysis failed (orbitals)')
#end if
#end try
if len(orbitals)>0:
self.orbitals = orbitals
#end if
#end if
# try to get the mulliken/lowdin populations in each ao
if 'counts' in self:
ao_populations = obj()
self.read_ao_populations(log,ao_populations)
try:
self.read_ao_populations(log,ao_populations)
except:
if self.info.exit:
self.error('log file analysis failed (ao populations)')
#end if
#end try
if len(ao_populations)>0:
self.ao_populations = ao_populations
#end if
#end if
#end def analyze_log
def read_energy_components(self,log,energy):
if log!=None and log.seek('ENERGY COMPONENTS',0)!=-1:
for n in xrange(18):
line = log.readline()
if '=' in line and 'ENERGY' in line:
nameline,value = line.split('=')
tokens = nameline.lower().split()
name = ''
for token in tokens:
if token!='energy':
name += token.replace('-','_')+'_'
#end if
#end for
name = name[:-1]
value = float(value.strip())
energy[name]=value
#end if
#end for
#end if
#end def read_energy_components
def read_orbitals(self,log,orbs,spec,header):
success = True
cao_tot = self.counts.cao
mos_tot = self.counts.mos
mos_found = 0
eigenvalue = []
symmetry = []
coefficients = []
have_basis = False
element = []
spec_index = []
angular = []
if log.seek(header,0)!=-1:
imax = mos_tot
jmax = 100
i=0
while mos_found<mos_tot and i<imax:
j=0
norbs = -1
while j<jmax:
line = log.readline()
if line.strip().replace(' ','').isdigit():
found_orbs = True
norbs = len(line.split())
break
#end if
j+=1
#end while
if j==jmax:
success=False
if self.info.exit:
self.error('could not find start of orbitals for {0}\nnumber of orbitals read successfully: {1}\nnumber of orbitals not read: {2}'.format(header,mos_found,mos_tot-mos_found))
else:
break
#end if
#end if
eigenvalue.extend(log.readtokens())
symmetry.extend(log.readtokens())
coeff = []
for icao in xrange(cao_tot):
tokens = log.readtokens()
if not have_basis:
e = tokens[1]
element.append(e[0].upper()+e[1:].lower())
spec_index.append(tokens[2])
angular.append(tokens[3])
#end if
coeff.append(tokens[4:])
#end for
coefficients.extend(array(coeff,dtype=float).T)
if not have_basis:
stype = []
ptype = []
dtype = []
ftype = []
for ia in xrange(len(angular)):
a = angular[ia].lower()
if a in GamessAnalyzer.stypes:
stype.append(ia)
elif a in GamessAnalyzer.ptypes:
ptype.append(ia)
elif a in GamessAnalyzer.dtypes:
dtype.append(ia)
elif a in GamessAnalyzer.ftypes:
ftype.append(ia)
elif self.info.exit:
self.error('unrecognized angular type: {0}'.format(angular[ia]))
#end if
#end for
#end if
have_basis = True
mos_found = len(eigenvalue)
i+=1
#end while
if i==imax:
success=False
if self.info.exit:
self.error('orbital read failed for {0}\nnumber of orbitals read successfully: {1}\nnumber of orbitals not read: {2}'.format(header,mos_found,mos_tot-mos_found))
#end if
if success:
orbs[spec] = obj(
eigenvalue = array(eigenvalue ,dtype=float),
symmetry = array(symmetry ,dtype=str ),
coefficients = array(coefficients,dtype=float),
basis = obj(
element = array(element ,dtype=str),
spec_index = array(spec_index,dtype=int),
angular = array(angular ,dtype=str),
stype = array(stype ,dtype=int),
ptype = array(ptype ,dtype=int),
dtype = array(dtype ,dtype=int),
ftype = array(ftype ,dtype=int),
)
)
#end if
return success
else:
return False
#end if
#end def read_orbitals
def read_ao_populations(self,log,ao_populations):
cao_tot = self.counts.cao
if log.seek('-- POPULATIONS IN EACH AO --',0)!=-1:
log.readline()
log.readline()
mulliken = []
lowdin = []
element = []
spec_index = []
angular = []
linds = obj()
for l in GamessAnalyzer.lset_full:
linds[l]=[]
#end for
for icao in xrange(cao_tot):
tokens = log.readtokens()
e = tokens[1]
element.append(e[0].upper()+e[1:].lower())
if len(tokens)==6:
spec_index.append(tokens[2])
angular.append(tokens[3])
mulliken.append(tokens[4])
lowdin.append(tokens[5])
elif len(tokens)==5:
spec_index.append(tokens[2][:-4])
angular.append(tokens[2][-4:])
mulliken.append(tokens[3])
lowdin.append(tokens[4])
#end if
#end for
for ia in xrange(len(angular)):
a = angular[ia].lower()
if a in GamessAnalyzer.lxyz_reverse:
l = GamessAnalyzer.lxyz_reverse[a]
linds[l].append(ia)
elif self.info.exit:
self.error('unrecognized angular type: {0}'.format(angular[ia]))
#end if
#end for
mulliken = array(mulliken,dtype=float)
lowdin = array(lowdin ,dtype=float)
for l,lind in linds.iteritems():
linds[l] = array(lind,dtype=int)
#end for
mulliken_shell = []
lowdin_shell = []
shell = obj()
n=0
for l in GamessAnalyzer.lset_full:
if l in GamessAnalyzer.lxyz:
lind = linds[l]
nxyz = len(GamessAnalyzer.lxyz[l])
nshell = len(lind)/nxyz
shellinds = []
for ns in range(nshell):
inds = lind[ns*nxyz:(ns+1)*nxyz]
mulliken_shell.append(mulliken[inds].sum())
lowdin_shell.append(lowdin[inds].sum())
shellinds.append(n)
n+=1
#end for
shell[l] = array(shellinds,dtype=int)
#end if
#end for
mulliken_shell = array(mulliken_shell)
lowdin_shell = array(lowdin_shell)
mulliken_angular = obj()
lowdin_angular = obj()
for l,shellinds in shell.iteritems():
mulliken_angular[l] = mulliken_shell[shellinds].sum()
lowdin_angular[l] = lowdin_shell[shellinds].sum()
#end for
basis = obj(
element = array(element ,dtype=str),
spec_index = array(spec_index,dtype=int),
angular = array(angular ,dtype=str),
)
basis.transfer_from(linds)
ao_populations.set(
mulliken = mulliken,
lowdin = lowdin,
mulliken_shell = mulliken_shell,
lowdin_shell = lowdin_shell,
mulliken_angular = mulliken_angular,
lowdin_angular = lowdin_angular,
basis = basis,
shell = shell,
)
#end if
#end def read_ao_populations
def analyze_punch(self):
# read the punch file
try:

View File

@ -141,7 +141,7 @@ class obj(AllAbilities):
if header==None:
header = self.__class__.__name__
#end if
self.log(header+post_header)
self.log('\n'+header+post_header)
self.log(pad+message.replace('\n','\n'+pad))
if exit:
self.log(' exiting.\n')
@ -157,7 +157,7 @@ class obj(AllAbilities):
if header==None:
header=self.__class__.__name__
#end if
self.log(header+post_header)
self.log('\n'+header+post_header)
self.log(pad+message.replace('\n','\n'+pad))
#end def warn

View File

@ -585,12 +585,12 @@ def generate_physical_system(**kwargs):
#end if
#end for
type = kwargs['type']
if type=='atom' or type=='dimer':
if type=='atom' or type=='dimer' or type=='trimer':
del kwargs['kshift']
del kwargs['tiling']
if not 'units' in kwargs:
kwargs['units'] = 'B'
#end if
#if not 'units' in kwargs:
# kwargs['units'] = 'B'
##end if
tiling = None
else:
tiling = kwargs['tiling']

View File

@ -2978,6 +2978,20 @@ class Structure(Sobj):
lcur = 7
#end if
species = elem
# relabel species that have multiple occurances
sset = set(species)
for spec in sset:
if species.count(spec)>1:
cnt=0
for n in range(len(species)):
specn = species[n]
if specn==spec:
cnt+=1
species[n] = specn+str(cnt)
#end if
#end for
#end if
#end for
elem = []
for i in range(len(counts)):
elem.extend(counts[i]*[species[i]])
@ -4000,6 +4014,8 @@ def generate_structure(type='crystal',*args,**kwargs):
s = generate_atom_structure(*args,**kwargs)
elif type=='dimer':
s = generate_dimer_structure(*args,**kwargs)
elif type=='trimer':
s = generate_trimer_structure(*args,**kwargs)
elif type=='jellium':
s = generate_jellium_structure(*args,**kwargs)
elif type=='empty':
@ -4007,7 +4023,7 @@ def generate_structure(type='crystal',*args,**kwargs):
elif type=='basic':
s = Structure(*args,**kwargs)
else:
Structure.class_error(str(type)+' is not a valid structure type\n options are crystal, defect, or atom')
Structure.class_error(str(type)+' is not a valid structure type\n options are crystal, defect, atom, dimer, trimer, jellium, empty, or basic')
#end if
return s
#end def generate_structure
@ -4016,6 +4032,9 @@ def generate_structure(type='crystal',*args,**kwargs):
def generate_atom_structure(atom=None,units='A',Lbox=None,skew=0,axes=None,kgrid=(1,1,1),kshift=(0,0,0),struct_type=Structure):
if atom is None:
Structure.class_error('atom must be provided','generate_atom_structure')
#end if
if Lbox!=None:
axes = [[Lbox*(1-skew),0,0],[0,Lbox,0],[0,0,Lbox*(1+skew)]]
#end if
@ -4030,6 +4049,9 @@ def generate_atom_structure(atom=None,units='A',Lbox=None,skew=0,axes=None,kgrid
def generate_dimer_structure(dimer=None,units='A',separation=None,Lbox=None,skew=0,axes=None,kgrid=(1,1,1),kshift=(0,0,0),struct_type=Structure,axis='x'):
if dimer is None:
Structure.class_error('dimer atoms must be provided to construct dimer','generate_dimer_structure')
#end if
if separation is None:
Structure.class_error('separation must be provided to construct dimer','generate_dimer_structure')
#end if
@ -4055,6 +4077,69 @@ def generate_dimer_structure(dimer=None,units='A',separation=None,Lbox=None,skew
#end def generate_dimer_structure
def generate_trimer_structure(trimer=None,units='A',separation=None,angle=None,Lbox=None,skew=0,axes=None,kgrid=(1,1,1),kshift=(0,0,0),struct_type=Structure,axis='x',axis2='y',angular_units='degrees'):
if trimer is None:
Structure.class_error('trimer atoms must be provided to construct trimer','generate_trimer_structure')
#end if
if separation is None:
Structure.class_error('separation must be provided to construct trimer','generate_trimer_structure')
#end if
if len(separation)!=2:
Structure.class_error('two separation distances (atom1-atom2,atom1-atom3) must be provided to construct trimer\nyou provided {0} separation distances'.format(len(separation)),'generate_trimer_structure')
#end if
if angle is None:
Structure.class_error('angle must be provided to construct trimer','generate_trimer_structure')
#end if
if angular_units=='degrees':
angle *= pi/180
elif not angle.startswith('rad'):
Structure.class_error('angular units must be degrees or radians\nyou provided: {0}'.format(angle),'generate_trimer_structure')
#end if
if axis==axis2:
Structure.class_error('axis and axis2 must be different to define the trimer plane\nyou provided {0} for both'.format(axis),'generate_trimer_structure')
#end if
if Lbox!=None:
axes = [[Lbox*(1-skew),0,0],[0,Lbox,0],[0,0,Lbox*(1+skew)]]
#end if
p1 = [0,0,0]
if axis=='x':
p2 = [separation[0],0,0]
elif axis=='y':
p2 = [0,separation[0],0]
elif axis=='z':
p2 = [0,0,separation[0]]
else:
Structure.class_error('trimer bond1 (atom2-atom1) orientation axis must be x,y,z\n you provided: {0}'.format(axis),'generate_trimer_structure')
#end if
r = separation[1]
c = cos(angle)
s = sin(angle)
axpair = axis+axis2
if axpair=='xy':
p3 = [r*c,r*s,0]
elif axpair=='yx':
p3 = [r*s,r*c,0]
elif axpair=='yz':
p3 = [0,r*c,r*s]
elif axpair=='zy':
p3 = [0,r*s,r*c]
elif axpair=='zx':
p3 = [r*s,0,r*c]
elif axpair=='xz':
p3 = [r*c,0,r*s]
else:
Structure.class_error('trimer bond2 (atom3-atom1) orientation axis must be x,y,z\n you provided: {0}'.format(axis2),'generate_trimer_structure')
#end if
if axes is None:
s = Structure(elem=trimer,pos=[p1,p2,p3],units=units)
else:
s = Structure(elem=trimer,pos=[p1,p2,p3],axes=axes,kgrid=kgrid,kshift=kshift,units=units)
s.center_molecule()
#end if
return s
#end def generate_trimer_structure
def generate_jellium_structure(*args,**kwargs):
return Jellium(*args,**kwargs)
#end def generate_jellium_structure

View File

@ -40,10 +40,11 @@
import os
from numpy import array,abs,empty,ndarray
from generic import obj
from periodic_table import is_element
from simulation import SimulationInput
from structure import interpolate_structures,Structure
from physical_system import PhysicalSystem
from developer import DevBase
from developer import DevBase,error
from debug import *
lcs = ls
@ -1008,7 +1009,11 @@ class Poscar(VFormattedFile):
#end for
if self.elem!=None:
for e in self.elem:
text += e+' '
iselem,symbol = is_element(e,symbol=True)
if not iselem:
self.error('{0} is not an element'.format(e))
#end if
text += symbol+' '
#end for
text += '\n'
#end if
@ -1276,7 +1281,7 @@ class VaspInput(SimulationInput,Vobj):
#end def write
def incorporate_system(self,system,incorp_kpoints=True):
def incorporate_system(self,system,incorp_kpoints=True,coord='cartesian'):
structure = system.structure
# assign kpoints
@ -1300,8 +1305,15 @@ class VaspInput(SimulationInput,Vobj):
poscar.axes = s.axes
poscar.elem = species
poscar.elem_count = species_count
poscar.coord = 'cartesian'
poscar.pos = s.pos
if coord=='cartesian':
poscar.coord = 'cartesian'
poscar.pos = s.pos
elif coord=='direct':
poscar.coord = 'direct'
poscar.pos = s.pos_unit()
else:
self.error('coord must be either direct or cartesian\nyou provided: {0}'.format(coord))
#end if
if s.frozen!=None:
poscar.dynamic = s.frozen==False
#end if
@ -1327,10 +1339,13 @@ class VaspInput(SimulationInput,Vobj):
#end for
ordered_pseudos = []
for element in species:
if not element in pseudo_map:
self.error('pseudopotential for element {0} not found\nelements present: {1}\n'.format(element,sorted(pseudo_map.keys())))
iselem,symbol = is_element(element,symbol=True)
if not iselem:
self.error('{0} is not an element'.format(element))
elif not symbol in pseudo_map:
self.error('pseudopotential for element {0} not found\nelements present: {1}'.format(symbol,sorted(pseudo_map.keys())))
#end if
ordered_pseudos.append(pseudo_map[element])
ordered_pseudos.append(pseudo_map[symbol])
#end for
#end if
self.potcar = Potcar(VaspInput.pseudo_dir,ordered_pseudos)
@ -1453,7 +1468,8 @@ generate_any_defaults = obj(
system = None,
pseudos = None,
neb = None,
neb_args = obj()
neb_args = obj(),
coord = 'cartesian'
)
def generate_any_vasp_input(**kwargs):
@ -1500,7 +1516,7 @@ def generate_any_vasp_input(**kwargs):
# incorporate system information
species = None
if vf.system!=None:
species = vi.incorporate_system(vf.system,gen_kpoints)
species = vi.incorporate_system(vf.system,gen_kpoints,vf.coord)
#end if
# set potcar
@ -1557,7 +1573,7 @@ def generate_any_vasp_input(**kwargs):
def generate_poscar(structure):
def generate_poscar(structure,coord='cartesian'):
s = structure.copy()
s.change_units('A')
species,species_count = s.order_by_species()
@ -1566,8 +1582,15 @@ def generate_poscar(structure):
poscar.axes = s.axes
poscar.elem = species
poscar.elem_count = species_count
poscar.coord = 'cartesian'
poscar.pos = s.pos
if coord=='cartesian':
poscar.coord = 'cartesian'
poscar.pos = s.pos
elif coord=='direct':
poscar.coord = 'direct'
poscar.pos = s.pos_unit()
else:
error('coord must be either direct or cartesian\nyou provided: {0}'.format(coord),'generate_poscar')
#end if
if s.frozen!=None:
poscar.dynamic = s.frozen==False
#end if