-
Notifications
You must be signed in to change notification settings - Fork 0
/
exclusions.py
72 lines (52 loc) · 2.8 KB
/
exclusions.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
import simtk.openmm as mm
import simtk.openmm.app as app
import numpy as np
from .constants import ONE_4PI_EPS0
def generateExclusion(self, force, forceName, param):
"""Create the exclusions for the Force object. Definition of exclusions will be inherited
from the nonbondedForce object."""
for (idx1, idx2, q_excl, sig_excl, eps_excl) in param['nbexclparam']:
force.addExclusion(idx1, idx2)
return force
def checkExclusion(self, listPair, param):
"""Check the exclusions in nonbondedForce object are overlapped with the atom in listPair"""
listexcl = []
for (idx1, idx2, q_excl, sig_excl, eps_excl) in param['nbexclparam']:
if ((idx1 in listPair[0]) and (idx2 in listPair[1])) or ((idx1 in listPair[1]) and (idx2 in listPair[0])):
listexcl.append([idx1,idx2,q_excl,sig_excl,eps_excl])
return listexcl
def createExclusionBondForce(self, listPair, forceName, method, param, tag=''):
"""Create the exclusion bonds for the CustomBondForce object. Usually 1-4 bond is selected."""
listexcl = self.checkExclusion(listPair, param)
if not listexcl:
return None
else:
print('Exclusions for '+forceName+' are now considered.')
if tag == 'PME':
exclexpr = "ONE_4PI_EPS0*chargeprod_excl/r - ONE_4PI_EPS0*chargeprod_nat*erf(alpha_ewald*r)/r;"
exclexpr += "ONE_4PI_EPS0 = {:f};".format(ONE_4PI_EPS0)
exclexpr += "alpha_ewald = {:f};".format(method['alpha_ewald'])
exclforce = mm.CustomBondForce(exclexpr)
exclforce.addPerBondParameter('chargeprod_excl')
exclforce.addPerBondParameter('chargeprod_nat')
for (idx1, idx2, q_excl, sig_excl, eps_excl) in param['nbexclparam']:
q1 = param['elecparam'][idx1]
q2 = param['elecparam'][idx2]
exclforce.addBond(idx1, idx2, [q_excl, q1*q2])
elif tag in ['NoCutoff','Cutoff']:
exclexpr = "ONE_4PI_EPS0*chargeprod_excl/r;"
exclexpr += "ONE_4PI_EPS0 = {:f};".format(ONE_4PI_EPS0)
exclforce = mm.CustomBondForce(exclexpr)
exclforce.addPerBondParameter('chargeprod_excl')
for (idx1, idx2, q_excl, sig_excl, eps_excl) in param['nbexclparam']:
exclforce.addBond(idx1, idx2, [q_excl])
elif tag == 'LennardJones':
exclexpr = "4*epsilon*((sigma/r)^12 - (sigma/r)^6);"
exclforce = mm.CustomBondForce(exclexpr)
exclforce.addPerBondParameter('sigma')
exclforce.addPerBondParameter('epsilon')
for (idx1, idx2, q_excl, sig_excl, eps_excl) in param['ljexclparam']:
exclforce.addBond(idx1, idx2, [sig_excl, eps_excl])
self.groupForce.append(exclforce)
self.groupForceName.append(forceName+'_excl')
return exclforce