Commit 03e33772 authored by ChayaSt's avatar ChayaSt
Browse files

using sigma instead of rmin

parent 73183c61
...@@ -1759,10 +1759,10 @@ class LennardJonesGenerator(object): ...@@ -1759,10 +1759,10 @@ class LennardJonesGenerator(object):
if None not in types: if None not in types:
type1 = types[0][0] type1 = types[0][0]
type2 = types[1][0] type2 = types[1][0]
emin = _convertParameterToNumber(parameters['emin']) epsilon = _convertParameterToNumber(parameters['epsilon'])
rmin = _convertParameterToNumber(parameters['rmin']) sigma = _convertParameterToNumber(parameters['sigma'])
self.nbfix_types[(type1, type2)] = [emin, rmin] self.nbfix_types[(type1, type2)] = [sigma, epsilon]
self.nbfix_types[(type2, type1)] = [emin, rmin] self.nbfix_types[(type2, type1)] = [sigma, epsilon]
def registerLennardJones(self, parameters): def registerLennardJones(self, parameters):
self.lj_types.registerAtom(parameters) self.lj_types.registerAtom(parameters)
...@@ -1817,20 +1817,20 @@ class LennardJonesGenerator(object): ...@@ -1817,20 +1817,20 @@ class LennardJonesGenerator(object):
for n in range(num_lj_types): for n in range(num_lj_types):
pair = (lj_type_list[m], lj_type_list[n]) pair = (lj_type_list[m], lj_type_list[n])
if pair in self.nbfix_types: if pair in self.nbfix_types:
wdij = self.nbfix_types[pair][0] wdij = self.nbfix_types[pair][1]
rij = self.nbfix_types[pair][1] rij = self.nbfix_types[pair][0] * 2**(1.0/6) / 2
rij6 = rij**6 rij6 = rij**6
acoef[m+num_lj_types*n] = math.sqrt(wdij) * rij6 acoef[m+num_lj_types*n] = wdij * rij6**2
bcoef[m+num_lj_types*n] = 2 * wdij * rij6 bcoef[m+num_lj_types*n] = 2 * wdij * rij6
continue continue
else: else:
rij = (reverse_map[m][0] + reverse_map[n][0]) * 2**(1.0/6) / 2 rij = (reverse_map[m][0] + reverse_map[n][0]) * 2**(1.0/6) / 2
rij6 = rij**6 rij6 = rij**6
wdij = math.sqrt(reverse_map[m][-1] * reverse_map[n][-1]) wdij = math.sqrt(reverse_map[m][-1] * reverse_map[n][-1])
acoef[m+num_lj_types*n] = math.sqrt(wdij) * rij6 acoef[m+num_lj_types*n] = wdij * rij6**2
bcoef[m+num_lj_types*n] = 2 * wdij * rij6 bcoef[m+num_lj_types*n] = 2 * wdij * rij6
self.force = mm.CustomNonbondedForce('a^2/r^12-b/r6; r6=r2*r2*r2; r2=r^2; ' self.force = mm.CustomNonbondedForce('a/r^12-b/r^6;'
'a=acoef(type1, type2);' 'a=acoef(type1, type2);'
'b=bcoef(type1, type2)') 'b=bcoef(type1, type2)')
self.force.addTabulatedFunction('acoef', self.force.addTabulatedFunction('acoef',
......
...@@ -759,7 +759,7 @@ class AmoebaTestForceField(unittest.TestCase): ...@@ -759,7 +759,7 @@ class AmoebaTestForceField(unittest.TestCase):
<LennardJonesForce lj14scale="1.0"> <LennardJonesForce lj14scale="1.0">
<Atom type="CLA" sigma="0.404468018036" epsilon="0.6276"/> <Atom type="CLA" sigma="0.404468018036" epsilon="0.6276"/>
<Atom type="SOD" sigma="0.251367073323" epsilon="0.1962296"/> <Atom type="SOD" sigma="0.251367073323" epsilon="0.1962296"/>
<NBFixPair type1="CLA" type2="SOD" emin="0.350933" rmin="0.3731"/> <NBFixPair type1="CLA" type2="SOD" sigma="0.664788623476" epsilon="0.350933"/>
</LennardJonesForce> </LennardJonesForce>
</ForceField> """ </ForceField> """
ff = ForceField(StringIO(xml)) ff = ForceField(StringIO(xml))
......
Markdown is supported
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