Commit 1f01ce9e authored by Berk Onat's avatar Berk Onat
Browse files

Correction for listing at interaction_type and atom_type definitions.

parent 4a891d98
......@@ -163,6 +163,28 @@ class MapDictionary(dict):
return [val.metaName for val in self.__dict__.values()]
def get_unitDict(keyname):
""" Unit dictionary for convertions
Unit names will be converted to values.
When defining units in translator dictionary,
the unit names in the dictionary should be used.
The unit convertion values are written for SI units.
If you would like to change it, just add another key
to the dictionary and change the key at parser base.
Usage:
Natively support python language in definitions.
You can use any python math operator and function.
Moreover, you can use space or - for multiplication
and ^ for power of values.
Example:
kilogram/meter^2 can be written as
kilo-gram/meter^2 or kilo-gram/meter**2
and will be calculated as
kilo*gram/meter**2
For a quick check for AMBER units please see:
http://ambermd.org/Questions/units.html
"""
unitDict = {
"si" : {
"meter" : "1.0",
......@@ -184,11 +206,10 @@ def get_unitDict(keyname):
"atto" : "1.0e-18",
"erg" : "1.0e-7",
"dyne" : "1.0e-5",
"barye" : "1.0e-1",
"bar" : "1.0e-1",
"angstrom" : "1.0e-10",
"kcal" : "4184.096739614824",
"mole" : "0.602213737699784e24",
"mol" : "0.602213737699784e24",
"atmosphere" : "1.01325e5",
"electron" : "1.602176565e-19",
"atomicmassunit" : "1.66054e-27",
......@@ -712,6 +733,7 @@ def get_updateDictionary(self, defname):
'value' : 'pres0'}
],
lookupdict=self.cntrlDict,
valtype='float',
#autoSections=True,
activeSections=['settings_barostat']
),
......@@ -722,6 +744,7 @@ def get_updateDictionary(self, defname):
'value' : 'taup'}
],
lookupdict=self.cntrlDict,
valtype='float',
#autoSections=True,
activeSections=['settings_barostat']
),
......@@ -742,12 +765,14 @@ def get_updateDictionary(self, defname):
'value' : 'dt'}
],
lookupdict=self.cntrlDict,
valtype='float',
#autoSections=True,
activeSections=['settings_integrator']
),
'x_amber_number_of_steps_requested' : MetaInfoMap(startpage,
depends=[{'value' : 'nstlim'}],
lookupdict=self.cntrlDict,
valtype='int',
#autoSections=True,
activeSections=['settings_integrator']
),
......@@ -774,6 +799,7 @@ def get_updateDictionary(self, defname):
'value' : 'temp0'}
],
lookupdict=self.cntrlDict,
valtype='float',
#autoSections=True,
activeSections=['settings_thermostat']
),
......@@ -789,6 +815,7 @@ def get_updateDictionary(self, defname):
'value' : 'gamma_ln'},
],
lookupdict=self.cntrlDict,
valtype='float',
#autoSections=True,
activeSections=['settings_thermostat']
),
......@@ -798,6 +825,7 @@ def get_updateDictionary(self, defname):
'value' : 'gamma_ln'}
],
lookupdict=self.cntrlDict,
valtype='float',
#autoSections=True,
activeSections=['settings_thermostat']
),
......@@ -830,23 +858,25 @@ def get_updateDictionary(self, defname):
'energy_current' : MetaInfoMap(startpage),
'energy_electrostatic' : MetaInfoMap(startpage,
depends=[{'value' : 'EELEC'}],
valtype='float',
unitdict=self.unitDict,
unit='electron-volt',
unit='kcal/mol',
lookupdict=self.mddataDict
),
'energy_free_per_atom' : MetaInfoMap(startpage),
'energy_free' : MetaInfoMap(startpage),
'energy_method_current' : MetaInfoMap(startpage,
depends=[{'assign' : 'Force Field'}],
lookupdict=self.mddataDict
),
#'energy_method_current' : MetaInfoMap(startpage,
# depends=[{'assign' : 'Force Field'}],
# lookupdict=self.mddataDict
# ),
'energy_T0_per_atom' : MetaInfoMap(startpage),
'energy_total_T0_per_atom' : MetaInfoMap(startpage),
'energy_total_T0' : MetaInfoMap(startpage),
'energy_total' : MetaInfoMap(startpage,
depends=[{'value' : 'Etot'}],
valtype='float',
unitdict=self.unitDict,
unit='electron-volt',
unit='kcal/mol',
lookupdict=self.mddataDict
),
'hessian_matrix' : MetaInfoMap(startpage),
......@@ -869,8 +899,9 @@ def get_updateDictionary(self, defname):
'energy_van_der_Waals_value' : MetaInfoMap(startpage,
depends=[{'value' : 'VDWAALS'}],
lookupdict=self.mddataDict,
valtype='float',
unitdict=self.unitDict,
unit='electron-volt',
unit='kcal/mol',
#autoSections=True,
activeSections=['section_energy_van_der_Waals']
),
......@@ -890,7 +921,7 @@ def get_updateDictionary(self, defname):
depends=[{'store' : 'RESTRAINT'}],
valtype='float',
unitdict=self.unitDict,
unit='electron-volt',
unit='kcal/mol',
lookupdict=self.mddataDict
),
'frame_sequence_continuation_kind' : MetaInfoMap(startpage),
......@@ -905,7 +936,7 @@ def get_updateDictionary(self, defname):
depends=[{'store' : 'EKtot'}],
valtype='float',
unitdict=self.unitDict,
unit='electron-volt',
unit='kcal/mol',
lookupdict=self.mddataDict
),
'frame_sequence_local_frames_ref' : MetaInfoMap(startpage),
......@@ -919,7 +950,7 @@ def get_updateDictionary(self, defname):
depends=[{'store' : 'EPtot'}],
valtype='float',
unitdict=self.unitDict,
unit='electron-volt',
unit='kcal/mol',
lookupdict=self.mddataDict
),
'frame_sequence_pressure_frames' : MetaInfoMap(startpage,
......@@ -1037,6 +1068,7 @@ def get_updateDictionary(self, defname):
'local_rotations' : MetaInfoMap(startpage),
'number_of_atoms' : MetaInfoMap(startpage,
depends=[{'value' : 'NATOM'}],
valtype='int',
lookupdict=self.parmDict
),
'number_of_sites' : MetaInfoMap(startpage),
......@@ -1067,6 +1099,7 @@ def get_updateDictionary(self, defname):
configuration_core = {
'number_of_electrons' : MetaInfoMap(startpage,
value=0,
valtype='float',
),
'atom_labels' : MetaInfoMap(startpage,
#subfunction=self.system_atom_labels()
......@@ -1116,13 +1149,15 @@ def get_updateDictionary(self, defname):
'molecule_to_molecule_type_map' : MetaInfoMap(startpage),
'number_of_topology_atoms' : MetaInfoMap(startpage,
depends=[{'value' : 'NATOM'}],
valtype='int',
lookupdict=self.parmDict
),
'number_of_topology_molecules' : MetaInfoMap(startpage,
subfunction=self.topology_num_topo_mol
subfunction=self.topology_num_topo_mol,
valtype='int',
),
'topology_force_field_name' : MetaInfoMap(startpage,
value='Amber Force Field',
value='AMBER Force Field',
)
}
......
......@@ -56,8 +56,8 @@ class AMBERParser(AmberC.AMBERParserBase):
'x_amber_geometry_optimization_cdetect': CachingLevel.Cache,
'x_amber_mdin_finline': CachingLevel.Ignore,
'x_amber_mdin_wt': CachingLevel.Ignore,
'x_amber_section_input_output_files': CachingLevel.Ignore,
'x_amber_single_configuration_calculation_detect': CachingLevel.Cache,
#'x_amber_single_configuration_calculation': CachingLevel.Cache,
}
for name in self.metaInfoEnv.infoKinds:
metaInfo = self.metaInfoEnv.infoKinds[name]
......@@ -82,6 +82,9 @@ class AMBERParser(AmberC.AMBERParserBase):
or name.startswith('section_single_configuration_calculation')
):
self.cachingLevelForMetaName[name] = CachingLevel.Cache
if name in self.extraDict.keys():
self.cachingLevelForMetaName[name] = CachingLevel.Ignore
def initialize_values(self):
"""Initializes the values of certain variables.
......@@ -230,7 +233,7 @@ class AMBERParser(AmberC.AMBERParserBase):
return False, None, itemdict
def topology_system_name(self, itemdict):
""" Function to generate data for atom_to_molecule
""" Function to generate data for system_name
"""
system_name = self.fName.split('.')[-1]
if self.topology:
......@@ -244,15 +247,17 @@ class AMBERParser(AmberC.AMBERParserBase):
residueList = self.topologyDict["resSeq"]
atomIndex = np.arange(len(residueList))
atom_to_residue = np.zeros((len(residueList), 2), dtype=int)
atom_to_residue[:,0] = atomIndex
atom_to_residue[:,1] = residueList
atom_to_residue[:,0] = atomIndex+1
atom_to_residue[:,1] = np.array(residueList)+1
return False, atom_to_residue, itemdict
else:
return False, None, itemdict
def topology_atom_type_and_interactions(self):
def topology_atom_type_and_interactions(self, backend, gIndex):
""" Function to generate data for atom_to_molecule
"""
sO = open_section
supB = backend.superBackend
if self.topology:
atom_type_list=list(set(self.topologyDict["name"]))
atom_type_dict = {}
......@@ -261,7 +266,7 @@ class AMBERParser(AmberC.AMBERParserBase):
radiusDict = {}
for ielm in range(len(atom_type_list)):
elm = atom_type_list[ielm]
atom_type_dict.update({elm : ielm})
atom_type_dict.update({elm : ielm+1})
for atom in self.topology.atoms:
if elm == atom.name:
massesDict.update({atom.name : atom.element.mass})
......@@ -272,62 +277,67 @@ class AMBERParser(AmberC.AMBERParserBase):
elementList = list(elementDict.values())
radiusList = list(radiusDict.values())
interNum = 0
interDict = {}
interTypeDict = {}
for ielm in range(len(atom_type_list)):
elm = atom_type_list[ielm]
bondList = []
typeList = []
bondid = 0
for bond in self.topology.bonds:
molname1, molatom1 = str(bond[0]).split('-')
molname2, molatom2 = str(bond[1]).split('-')
if (elm == str(molatom1) or elm == str(molatom2)):
bondList.append(list(self.topologyBonds[bondid]))
typeList.append(list([
atom_type_dict[str(molatom1)],
atom_type_dict[str(molatom2)]
]))
#str(molatom1)+':'+str(ielm+1),
#str(molatom2)+':'+str(ielm+1)
#]))
interDict.update({elm : bondList})
interTypeDict.update({elm : typeList})
bondid += 1
for ielm in range(len(atom_type_list)-1):
for jelm in range(ielm+1, len(atom_type_list)):
aelm = atom_type_list[ielm]
belm = atom_type_list[jelm]
bondList = []
typeList = []
bondid = 0
noninter = True
for bond in self.topology.bonds:
molname1, molatom1 = str(bond[0]).split('-')
molname2, molatom2 = str(bond[1]).split('-')
if((aelm == str(molatom1) and belm == str(molatom2)) or
(aelm == str(molatom2) and belm == str(molatom1))):
bondList.append(list(self.topologyBonds[bondid]))
interDict.update({interNum : bondList})
if noninter:
noninter = False
typeList.extend(list([
atom_type_dict[aelm],
atom_type_dict[belm]
]))
interTypeDict.update({interNum : typeList})
interNum += 1
bondid += 1
for elm in range(len(atom_type_list)):
with self.secOpen(self.superP, 'section_atom_type'):
with sO(supB, 'section_atom_type'):
### !!! ----------------------------------------------- !!!
### !!! Need to check the definition of atom type name.
### !!! Which one is better? C1, CA, C:1 or C !!!
### !!! ----------------------------------------------- !!!
# Atom name? Atom element?
self.superP.addValue('atom_type_name', str(atom_type_list[elm]) + ':' + str(elm+1))
supB.addValue('atom_type_name', str(atom_type_list[elm]) + ':' + str(elm+1))
# Atomic mass
self.superP.addValue('atom_type_mass', massesList[elm])
supB.addValue('atom_type_mass', massesList[elm])
# Atomic element/symbol
self.superP.addValue('x_amber_atom_type_element', elementList[elm])
supB.addValue('x_amber_atom_type_element', elementList[elm])
# Atomic van der Waals radius
self.superP.addValue('x_amber_atom_type_radius', radiusList[elm])
supB.addValue('x_amber_atom_type_radius', radiusList[elm])
# Atomic charge
#self.superP.addValue('atom_type_charge', atom_charge)
pass
for elm in atom_type_list:
with self.secOpen(self.superP, 'section_interaction'):
for inum in range(interNum):
with sO(supB, 'section_interaction'):
# atom indexes of bound pairs for a specific atom type
self.superP.addArrayValues('interaction_atoms', np.asarray(interDict[elm]))
supB.addArrayValues('interaction_atoms', np.asarray(interDict[inum]))
# number of bonds for this type
self.superP.addValue('number_of_interactions', len(interDict[elm]))
supB.addValue('number_of_interactions', len(interDict[inum]))
# number of atoms involved (2 for covalent bonds)
self.superP.addValue('number_of_atoms_per_interaction', len(interDict[elm][0]))
supB.addValue('number_of_atoms_per_interaction', len(interDict[inum][0]))
#if bondFunctional:
# self.superP.addValue('interaction_kind', bondFunctional) # functional form of the interaction
# this points to the relative section_atom_type
self.superP.addArrayValues('x_amber_interaction_atom_to_atom_type_ref',
np.asarray(interTypeDict[elm]))
supB.addArrayValues('x_amber_interaction_atom_to_atom_type_ref',
np.asarray(interTypeDict[inum]))
# interaction parameters for the functional
#self.superP.addValue('interaction_parameters', bondParameters)
......@@ -385,25 +395,7 @@ class AMBERParser(AmberC.AMBERParserBase):
The scala part will check the validity nevertheless. (Good to know :)
"""
# write trajectory
valuesDict = section.simpleValues
location = 'verbatim writeout of mdin',
# write settings of aims output from the parsed control.in
for k,v in section.simpleValues.items():
if (k.startswith('x_amber_mdin_') or
k.startswith('x_amber_mdout_') or
k.startswith('x_amber_parm_')):
backend.superBackend.addValue(k, v[-1])
# reset all variables
self.initialize_values()
# if self.MD:
# sampling_method = "molecular_dynamics"
## elif len(self.singleConfCalcs) > 1:
# pass # to do
# else:
# return
# check for geometry optimization convergence
#valuesDict = section.simpleValues
# samplingGIndex = backend.openSection("section_sampling_method")
# backend.addValue("sampling_method", sampling_method)
......@@ -416,6 +408,9 @@ class AMBERParser(AmberC.AMBERParserBase):
backend.addArrayValues("frame_sequence_local_frames_ref", np.asarray(self.singleConfCalcs))
backend.closeSection("section_frame_sequence", frameSequenceGIndex)
# reset all variables
self.initialize_values()
def onClose_x_amber_section_input_output_files(self, backend, gIndex, section):
"""Trigger called when x_amber_section_input_output_files is closed.
......@@ -510,7 +505,7 @@ class AMBERParser(AmberC.AMBERParserBase):
self.metaStorage.updateBackend(backend.superBackend,
startsection=['section_topology'],
autoopenclose=False)
self.topology_atom_type_and_interactions()
self.topology_atom_type_and_interactions(backend, gIndex)
if (gIndex is None or gIndex == -1 or gIndex == "-1"):
backend.superBackend.closeSection("section_topology", self.secTopologyGIndex)
......
......@@ -261,28 +261,65 @@ class Container(object):
postfunc = item.postfunction
storeValue, updateValue, item = postfunc(item)
if "valtype" in item:
if(isinstance(updateValue, list) or isinstance(updateValue, tuple)):
newUpdateValue = [eval(
item["valtype"]+"("+str(ival)+")"
) for ival in updateValue]
elif isinstance(updateValue, str):
newUpdateValue = eval(item["valtype"]+"("+str(updateValue)+")")
else:
newUpdateValue = updateValue
if updateValue is not None:
updateValue = self.convertToNumber(updateValue, item["valtype"])
if("unit" in item and "unitdict" in item):
if(isinstance(updateValue, list) or isinstance(updateValue, tuple)):
updateValue = [self.unitConverter(
ival, item["unit"], item["unitdict"]
) for ival in updateValue]
elif(isinstance(updateValue, str) and float(updateValue)):
updateValue = float(updateValue)
updateValue = self.unitConverter(
updateValue, item["unit"], item["unitdict"])
else:
updateValue = self.unitConverter(
updateValue, item["unit"], item["unitdict"])
if updateValue is not None:
updateValue = self.convertUnits(updateValue, item["unit"], item["unitdict"])
return storeValue, updateValue, localdict
def convertToNumber(self, updateValue, valtype):
if(isinstance(updateValue, list) or isinstance(updateValue, tuple)):
newUpdateValue = [eval(
valtype+"("+str(ival)+")"
) for ival in updateValue]
elif isinstance(updateValue, np.ndarray):
newUpdateValue = np.asarray([eval(
valtype+"("+str(ival)+")"
) for ival in updateValue])
elif self.is_number(updateValue):
newUpdateValue = eval(valtype+"("+str(updateValue)+")")
else:
# I hope you know what you are doing
newUpdateValue = updateValue
return newUpdateValue
def convertUnits(self, updateValue, unit, unitdict):
if(isinstance(updateValue, list) or isinstance(updateValue, tuple)):
updateValue = [self.unitConverter(
ival, unit, unitdict
) for ival in updateValue]
elif isinstance(updateValue, np.ndarray):
updateValue = np.asarray([self.unitConverter(
ival, unit, unitdict
) for ival in updateValue])
elif self.is_number(updateValue):
updateValue = self.convertToNumber(updateValue, "float")
updateValue = self.unitConverter(
updateValue, unit, unitdict)
else:
# I hope you know what you are doing
updateValue = self.unitConverter(
updateValue, unit, unitdict)
def unitConverter(self, updateValue, unit, unitdict):
""" Unit converter using definitions of units explicitly
The unit names are converted to numbers and the resulting
expression will be evaluated by python.
Ex.: unit = 'electron-volt/Angstrom^2'
will be converted to
unit = '1.602176565e-19*1.0/1.0e-10**2'
factor = eval(unit) = 16.02176565 Joule/meter^2
in SI units and the result will be calculated as follows:
output_value = input_value * factor
"""
newunit = unit.lower()
newunit = newunit.replace('-','*').replace(' ', '*').replace('^', "**")
for key,value in unitdict.items():
newunit = newunit.replace(str(key), str(value))
return float(updateValue) * float(eval(newunit))
def checkTestsDicts(self, item, localdict):
for depdict in item["depends"]:
deptests = depdict["test"]
......@@ -342,24 +379,6 @@ class Container(object):
return checkval, localdict
return None, localdict
def unitConverter(self, updateValue, unit, unitdict):
""" Unit converter using definitions of units explicitly
The unit names are converted to numbers and the resulting
expression will be evaluated by python.
Ex.: unit = 'electron-volt/Angstrom^2'
will be converted to
unit = '1.602176565e-19*1.0/1.0e-10**2'
factor = eval(unit) = 16.02176565 Joule/meter^2
in SI units and the result will be calculated as follows:
output_value = input_value * factor
"""
newunit = unit.lower()
newunit = newunit.replace('-','*').replace(' ', '*').replace('^', "**")
for key,value in unitdict.items():
newunit = newunit.replace(str(key), str(value))
return float(updateValue) * float(eval(newunit))
def findNameInLookupDict(self, metaname, lookupdict):
for item in lookupdict:
itemMap = lookupdict[item]
......@@ -430,6 +449,13 @@ class Container(object):
newValue = itemv["unitconverter"](self, itemv)
self.Storage.__dict__[itemk["val"]] = newvalue
def is_number(self, val):
try:
float(val)
return True
except ValueError:
return False
def __str__(self, caller=None, decorate='', color=None, printactive=None, onlynames=None):
string = ''
if onlynames is None:
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment