forcefield.py 164 KB
Newer Older
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
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
"""
forcefield.py: Constructs OpenMM System objects based on a Topology and an XML force field description
"""
__author__ = "Peter Eastman"
__version__ = "1.0"

import os
import itertools
import xml.etree.ElementTree as etree
from math import sqrt, cos
import simtk.openmm as mm
import simtk.unit as unit
import element as elem
from simtk.openmm.app import Topology

# Enumerated values for nonbonded method

NoCutoff = object()
CutoffNonPeriodic = object()
CutoffPeriodic = object()
Ewald = object()
PME = object()

# Enumerated values for constraint type

HBonds = object()
AllBonds = object()
HAngles = object()

# A map of functions to parse elements of the XML file.

parsers = {}

class ForceField(object):
    """A ForceField constructs OpenMM System objects based on a Topology."""

    def __init__(self, *files):
        """Load one or more XML files and create a ForceField object based on them.
        
        Parameters:
         - A list of XML files defining the force field.  Each entry may be an absolute file path, a path relative to the
           current working directory, or a path relative to this module's data subdirectory (for built in force fields).
        """
        self._atomTypes = {}
        self._templates = {}
        self._templateSignatures = {}
        self._atomClasses = {}
        self._forces = []
        for file in files:
            try:
                tree = etree.parse(file)
            except IOError:
                tree = etree.parse(os.path.join(os.path.dirname(__file__), 'data', file))
            root = tree.getroot()
            
            # Load the atom types.
            
            if tree.getroot().find('AtomTypes') is not None:
                for type in tree.getroot().find('AtomTypes').findall('Type'):
                    self._atomTypes[type.attrib['name']] = (type.attrib['class'], float(type.attrib['mass']), elem.get_by_symbol(type.attrib['element']))
            
            # Load the residue templates.
            
            if tree.getroot().find('Residues') is not None:
                for residue in root.find('Residues').findall('Residue'):
                    resName = residue.attrib['name']
                    template = ForceField._TemplateData(resName)
                    self._templates[resName] = template
                    for atom in residue.findall('Atom'):
                        template.atoms.append(ForceField._TemplateAtomData(atom.attrib['name'], atom.attrib['type'], self._atomTypes[atom.attrib['type']][2]))
                    for bond in residue.findall('Bond'):
                        b = (int(bond.attrib['from']), int(bond.attrib['to']))
                        template.bonds.append(b)
                        template.atoms[b[0]].bondedTo.append(b[1])
                        template.atoms[b[1]].bondedTo.append(b[0])
                    for bond in residue.findall('ExternalBond'):
                        b = int(bond.attrib['from'])
                        template.externalBonds.append(b)
                        template.atoms[b].externalBonds += 1
            for template in self._templates.values():
                template.signature = _createResidueSignature([atom.element for atom in template.atoms])
                sigString = _signatureToString(template.signature)
                if sigString in self._templateSignatures:
                    self._templateSignatures[sigString].append(template)
                else:
                    self._templateSignatures[sigString] = [template]
            
            # Build sets of every atom type belonging to each class
            
            for type in self._atomTypes:
                atomClass = self._atomTypes[type][0]
                if atomClass in self._atomClasses:
                    typeSet = self._atomClasses[atomClass]
                else:
                    typeSet = set()
                    self._atomClasses[atomClass] = typeSet
                typeSet.add(type)
            self._atomClasses[''] = self._atomTypes.keys()
            
            # Load force definitions
            
            for child in root:
                if child.tag in parsers:
                    parsers[child.tag](child, self)

    def _findAtomTypes(self, node, num):
        """Parse the attributes on an XML tag to find the set of atom types for each atom it involves."""
        types = []
        attrib = node.attrib
        for i in range(num):
            if num == 1:
                suffix = ''
            else:
                suffix = str(i+1)
            classAttrib = 'class'+suffix
            if classAttrib in attrib:
                if attrib[classAttrib] not in self._atomClasses:
                    return None # Unknown atom class
                types.append(self._atomClasses[attrib[classAttrib]])
            else:
                typeAttrib = 'type'+suffix
                if typeAttrib not in attrib or attrib[typeAttrib] not in self._atomTypes:
                    return None # Unknown atom type
                types.append([attrib[typeAttrib]])
        return types

    def _parseTorsion(self, node):
        """Parse the node defining a torsion."""
        types = self._findAtomTypes(node, 4)
        if types is None:
            return None
        torsion = PeriodicTorsion(types)
        attrib = node.attrib
        index = 1
        while 'phase%d'%index in attrib:
            torsion.periodicity.append(int(attrib['periodicity%d'%index]))
            torsion.phase.append(float(attrib['phase%d'%index]))
            torsion.k.append(float(attrib['k%d'%index]))
            index += 1
        return torsion
        
    class _SystemData:
        """Inner class used to encapsulate data about the system being created."""
        def __init__(self):
            self.atomType = {}
            self.atoms = []
            self.bonds = []
            self.angles = []
            self.propers = []
            self.impropers = []
            self.atomBonds = []
            self.isAngleConstrained = []

    class _TemplateData:
        """Inner class used to encapsulate data about a residue template definition."""
        def __init__(self, name):
            self.name = name
            self.atoms = []
            self.bonds = []
            self.externalBonds = []

    class _TemplateAtomData:
        """Inner class used to encapsulate data about an atom in a residue template definition."""
        def __init__(self, name, type, element):
            self.name = name
            self.type = type
            self.element = element
            self.bondedTo = []
            self.externalBonds = 0

    class _BondData:
        """Inner class used to encapsulate data about a bond."""
        def __init__(self, atom1, atom2):
            self.atom1 = atom1
            self.atom2 = atom2
            self.isConstrained = False
            self.length = 0.0

    def createSystem(self, topology, nonbondedMethod=NoCutoff, nonbondedCutoff=1.0*unit.nanometer,
                     constraints=None, rigidWater=True, **args):
        """Construct an OpenMM System representing a Topology with this force field.
        
        Parameters:
         - topology (Topology) The Topology for which to create a System
         - nonbondedMethod (object=NoCutoff) The method to use for nonbonded interactions.  Allowed values are
           NoCutoff, CutoffNonPeriodic, CutoffPeriodic, Ewald, or PME.
         - nonbondedCutoff (distance=1*nanometer) The cutoff distance to use for nonbonded interactions
         - constraints (object=None) Specifies which bonds angles should be implemented with constraints.
           Allowed values are None, HBonds, AllBonds, or HAngles.
         - rigidWater (boolean=True) If true, water molecules will be fully rigid regardless of the value passed for the constraints argument
         - Arbitrary additional keyword arguments may also be specified.  This allows extra parameters to be specified that are specific to
           particular force fields.
        Returns: the newly created System
        """
        
        # Record atom indices
        
        data = ForceField._SystemData()
        atomIndices = {}
        for index, atom in enumerate(topology.atoms()):
            data.atoms.append(atom)
            atomIndices[atom] = index

        # Make a list of all bonds
        
        for bond in topology.bonds():
            if bond[0] in atomIndices and bond[1] in atomIndices:
                data.bonds.append(ForceField._BondData(atomIndices[bond[0]], atomIndices[bond[1]]))

        # Record which atoms are bonded to each other atom
        
        bondedToAtom = []
        for i in range(len(data.atoms)):
            bondedToAtom.append(set())
            data.atomBonds.append([])
        for i in range(len(data.bonds)):
            bond = data.bonds[i]
            bondedToAtom[bond.atom1].add(bond.atom2)
            bondedToAtom[bond.atom2].add(bond.atom1)
            data.atomBonds[bond.atom1].append(i)
            data.atomBonds[bond.atom2].append(i)

        # Find the template matching each residue and assign atom types.
        
        for chain in topology.chains():
            for res in chain.residues():
                signature = _signatureToString(_createResidueSignature([atom.element for atom in res.atoms()]))
                template = None
                matches = None
                if signature in self._templateSignatures:
                    for t in self._templateSignatures[signature]:
                        matches = _matchResidue(res, t, bondedToAtom, atomIndices)
                        if matches is not None:
                            template = t
                            break
                if matches is None:
                    raise ValueError('No template found for residue %d (%s)' % (res.index+1, res.name))
                for atom, match in zip(res.atoms(), matches):
                    data.atomType[atom] = template.atoms[match].type

        # Create the System and add atoms
        
        sys = mm.System()
        for atom in topology.atoms():
            sys.addParticle(self._atomTypes[data.atomType[atom]][1])
        
        # Set periodic boundary conditions.
        
        boxSize = topology.getUnitCellDimensions()
        if boxSize is not None:
            sys.setDefaultPeriodicBoxVectors((boxSize[0], 0, 0), (0, boxSize[1], 0), (0, 0, boxSize[2]))
        elif nonbondedMethod is not NoCutoff and nonbondedMethod is not CutoffNonPeriodic:
            raise ValueError('Requested periodic boundary conditions for a Topology that does not specify periodic box dimensions')

        # Make a list of all unique angles
        
        uniqueAngles = set()
        for bond in data.bonds:
            for atom in bondedToAtom[bond.atom1]:
                if atom != bond.atom2:
                    if atom < bond.atom2:
                        uniqueAngles.add((atom, bond.atom1, bond.atom2))
                    else:
                        uniqueAngles.add((bond.atom2, bond.atom1, atom))
            for atom in bondedToAtom[bond.atom2]:
                if atom != bond.atom1:
                    if atom > bond.atom1:
                        uniqueAngles.add((bond.atom1, bond.atom2, atom))
                    else:
                        uniqueAngles.add((atom, bond.atom2, bond.atom1))
        data.angles = sorted(list(uniqueAngles))
        
        # Make a list of all unique proper torsions
        
        uniquePropers = set()
        for angle in data.angles:
            for atom in bondedToAtom[angle[0]]:
                if atom != angle[1]:
                    if atom < angle[2]:
                        uniquePropers.add((atom, angle[0], angle[1], angle[2]))
                    else:
                        uniquePropers.add((angle[2], angle[1], angle[0], atom))
            for atom in bondedToAtom[angle[2]]:
                if atom != angle[1]:
                    if atom > angle[0]:
                        uniquePropers.add((angle[0], angle[1], angle[2], atom))
                    else:
                        uniquePropers.add((atom, angle[2], angle[1], angle[0]))
        data.propers = sorted(list(uniquePropers))
        
        # Make a list of all unique improper torsions
        
        for atom in range(len(bondedToAtom)):
            bondedTo = bondedToAtom[atom]
            if len(bondedTo) > 2:
                for subset in itertools.combinations(bondedTo, 3):
                    data.impropers.append((atom, subset[0], subset[1], subset[2]))
        
        # Identify bonds that should be implemented with constraints
        
        if constraints == AllBonds or constraints == HAngles:
            for bond in data.bonds:
                bond.isConstrained = True
        elif constraints == HBonds:
            for bond in data.bonds:
                atom1 = data.atoms[bond.atom1]
                atom2 = data.atoms[bond.atom2]
                bond.isConstrained = atom1.name.startswith('H') or atom2.name.startswith('H')
        if rigidWater:
            for bond in data.bonds:
                atom1 = data.atoms[bond.atom1]
                atom2 = data.atoms[bond.atom2]
                if atom1.residue.name == 'HOH' and atom2.residue.name == 'HOH':
                    bond.isConstrained = True
        
        # Identify angles that should be implemented with constraints
        
        if constraints == HAngles:
            for angle in data.angles:
                atom1 = data.atoms[angle[0]]
                atom2 = data.atoms[angle[1]]
                atom3 = data.atoms[angle[2]]
                numH = 0
                if atom1.name.startswith('H'):
                    numH += 1
                if atom3.name.startswith('H'):
                    numH += 1
                data.isAngleConstrained.append(numH == 2 or (numH == 1 and atom2.name.startswith('O')))
        else:
            data.isAngleConstrained = len(data.angles)*[False]
        if rigidWater:
            for i in range(len(data.angles)):
                angle = data.angles[i]
                atom1 = data.atoms[angle[0]]
                atom2 = data.atoms[angle[1]]
                atom3 = data.atoms[angle[2]]
                if atom1.residue.name == 'HOH' and atom2.residue.name == 'HOH' and atom3.residue.name == 'HOH':
                    data.isAngleConstrained[i] = True

        # Add forces to the System
        
        for force in self._forces:
            force.createForce(sys, data, nonbondedMethod, nonbondedCutoff, args)
        return sys


def _createResidueSignature(elements):
    """Create a signature for a residue based on the elements of the atoms it contains."""
    counts = {}
    for element in elements:
        if element in counts:
            counts[element] += 1
        else:
            counts[element] = 1
    sig = []
    for c in counts:
        sig.append((c, counts[c]))
    sig.sort(key=lambda x: -x[0].mass)
    return sig


def _signatureToString(signature):
    """Convert the signature returned by _createResidueSignature() to a string."""
    s = ''
    for element, count in signature:
        s += element.symbol+str(count)
    return s


def _matchResidue(res, template, bondedToAtom, atomIndices):
    """Determine whether a residue matches a template and return a list of corresponding atoms.
    
    Parameters:
     - res (Residue) The residue to check
     - template (_TemplateData) The template to compare it to
     - bondedToAtom (list) Enumerates which other atoms each atom is bonded to
     - atomIndices (map) Maps from atoms to their indices in the System
    Returns: a list specifying which atom of the template each atom of the residue corresponds to,
    or None if it does not match the template
    """
    atoms = list(res.atoms())
    residueAtomBonds = []
    templateAtomBonds = []
    matches = len(atoms)*[0]
    hasMatch = len(atoms)*[False]
    
    # Translate from global to local atom indices, and record the bonds for each atom.
    
    renumberAtoms = {}
    for i in range(len(atoms)):
        renumberAtoms[atomIndices[atoms[i]]] = i
    bondedTo = []
    externalBonds = []
    for atom in atoms:
        bonds = [renumberAtoms[x] for x in bondedToAtom[atomIndices[atom]] if x in renumberAtoms]
        bondedTo.append(bonds)
        externalBonds.append(len([x for x in bondedToAtom[atomIndices[atom]] if x not in renumberAtoms]))
    if _findAtomMatches(atoms, template, bondedTo, externalBonds, matches, hasMatch, 0):
        return matches
    return None


def _findAtomMatches(atoms, template, bondedTo, externalBonds, matches, hasMatch, position):
    """This is called recursively from inside _matchResidue() to identify matching atoms."""
    if position == len(atoms):
        return True
    elem = atoms[position].element
    for i in range(len(atoms)):
        atom = template.atoms[i]
        if atom.element == elem and not hasMatch[i] and len(atom.bondedTo) == len(bondedTo[position]) and atom.externalBonds == externalBonds[position]:
            # See if the bonds for this identification are consistent
            
            allBondsMatch = all((bonded > position or matches[bonded] in atom.bondedTo for bonded in bondedTo[position]))
            if allBondsMatch:
                # This is a possible match, so trying matching the rest of the residue.
                
                matches[position] = i
                hasMatch[i] = True
                if _findAtomMatches(atoms, template, bondedTo, externalBonds, matches, hasMatch, position+1):
                    return True
                hasMatch[i] = False
    return False


# The following classes are generators that know how to create Force subclasses and add them to a System that is being
# created.  Each generator class must define two methods: 1) a static method that takes an etree Element and a ForceField,
# and returns the corresponding generator object; 2) a createForce() method that constructs the Force object and adds it
# to the System.  The static method should be added to the parsers map.


class HarmonicBondGenerator:
    """A HarmonicBondGenerator constructs a HarmonicBondForce."""
    
    def __init__(self):
        self.types1 = []
        self.types2 = []
        self.length = []
        self.k = []
    
    @staticmethod
    def parseElement(element, ff):
        generator = HarmonicBondGenerator()
        ff._forces.append(generator)
        for bond in element.findall('Bond'):
            types = ff._findAtomTypes(bond, 2)
            if types is not None:
                generator.types1.append(types[0])
                generator.types2.append(types[1])
                generator.length.append(float(bond.attrib['length']))
                generator.k.append(float(bond.attrib['k']))
    
    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
        existing = [sys.getForce(i) for i in range(sys.getNumForces())]
        existing = [f for f in existing if type(f) == mm.HarmonicBondForce]
        if len(existing) == 0:
            force = mm.HarmonicBondForce()
            sys.addForce(force)
        else:
            force = existing[0]
        for bond in data.bonds:
            type1 = data.atomType[data.atoms[bond.atom1]]
            type2 = data.atomType[data.atoms[bond.atom2]]
            for i in range(len(self.types1)):
                types1 = self.types1[i]
                types2 = self.types2[i]
                if (type1 in types1 and type2 in types2) or (type1 in types2 and type2 in types1):
                    bond.length = self.length[i]
                    if bond.isConstrained:
                        sys.addConstraint(bond.atom1, bond.atom2, self.length[i])
                    elif self.k[i] != 0:
                        force.addBond(bond.atom1, bond.atom2, self.length[i], self.k[i])
                    break

parsers["HarmonicBondForce"] = HarmonicBondGenerator.parseElement


class HarmonicAngleGenerator:
    """A HarmonicAngleGenerator constructs a HarmonicAngleForce."""
    
    def __init__(self):
        self.types1 = []
        self.types2 = []
        self.types3 = []
        self.angle = []
        self.k = []
    
    @staticmethod
    def parseElement(element, ff):
        generator = HarmonicAngleGenerator()
        ff._forces.append(generator)
        for angle in element.findall('Angle'):
            types = ff._findAtomTypes(angle, 3)
            if types is not None:
                generator.types1.append(types[0])
                generator.types2.append(types[1])
                generator.types3.append(types[2])
                generator.angle.append(float(angle.attrib['angle']))
                generator.k.append(float(angle.attrib['k']))
    
    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
        existing = [sys.getForce(i) for i in range(sys.getNumForces())]
        existing = [f for f in existing if type(f) == mm.HarmonicAngleForce]
        if len(existing) == 0:
            force = mm.HarmonicAngleForce()
            sys.addForce(force)
        else:
            force = existing[0]
        for (angle, isConstrained) in zip(data.angles, data.isAngleConstrained):
            type1 = data.atomType[data.atoms[angle[0]]]
            type2 = data.atomType[data.atoms[angle[1]]]
            type3 = data.atomType[data.atoms[angle[2]]]
            for i in range(len(self.types1)):
                types1 = self.types1[i]
                types2 = self.types2[i]
                types3 = self.types3[i]
                if (type1 in types1 and type2 in types2 and type3 in types3) or (type1 in types3 and type2 in types2 and type3 in types1):
                    if isConstrained:
                        # Find the two bonds that make this angle.
                        
                        bond1 = None
                        bond2 = None
                        for bond in data.atomBonds[angle[1]]:
                            atom1 = data.bonds[bond].atom1
                            atom2 = data.bonds[bond].atom2
                            if atom1 == angle[0] or atom2 == angle[0]:
                                bond1 = bond
                            elif atom1 == angle[2] or atom2 == angle[2]:
                                bond2 = bond
                        
                        # Compute the distance between atoms and add a constraint
                        
                        if bond1 is not None and bond2 is not None:
                            l1 = data.bonds[bond1].length
                            l2 = data.bonds[bond2].length
                            if l1 is not None and l2 is not None:
                                length = sqrt(l1*l1 + l2*l2 - 2*l1*l2*cos(self.angle[i]))
                                sys.addConstraint(angle[0], angle[2], length)
                    elif self.k[i] != 0:
                        force.addAngle(angle[0], angle[1], angle[2], self.angle[i], self.k[i])
                    break

parsers["HarmonicAngleForce"] = HarmonicAngleGenerator.parseElement


class PeriodicTorsion:
    """A PeriodicTorsion records the information for a periodic torsion definition."""

    def __init__(self, types):
        self.types1 = types[0]
        self.types2 = types[1]
        self.types3 = types[2]
        self.types4 = types[3]
        self.periodicity = []
        self.phase = []
        self.k = []

class PeriodicTorsionGenerator:
    """A PeriodicTorsionGenerator constructs a PeriodicTorsionForce."""
    
    def __init__(self):
        self.proper = []
        self.improper = []
    
    @staticmethod
    def parseElement(element, ff):
        generator = PeriodicTorsionGenerator()
        generator.ff = ff
        ff._forces.append(generator)
        for torsion in element.findall('Proper'):
            torsion = ff._parseTorsion(torsion)
            if torsion is not None:
                generator.proper.append(torsion)
        for torsion in element.findall('Improper'):
            torsion = ff._parseTorsion(torsion)
            if torsion is not None:
                generator.improper.append(torsion)
    
    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
        existing = [sys.getForce(i) for i in range(sys.getNumForces())]
        existing = [f for f in existing if type(f) == mm.PeriodicTorsionForce]
        if len(existing) == 0:
            force = mm.PeriodicTorsionForce()
            sys.addForce(force)
        else:
            force = existing[0]
        wildcard = self.ff._atomClasses['']
        for torsion in data.propers:
            type1 = data.atomType[data.atoms[torsion[0]]]
            type2 = data.atomType[data.atoms[torsion[1]]]
            type3 = data.atomType[data.atoms[torsion[2]]]
            type4 = data.atomType[data.atoms[torsion[3]]]
            match = None
            for tordef in self.proper:
                types1 = tordef.types1
                types2 = tordef.types2
                types3 = tordef.types3
                types4 = tordef.types4
                if (type2 in types2 and type3 in types3 and type4 in types4 and type1 in types1) or (type2 in types3 and type3 in types2 and type4 in types1 and type1 in types4):
                    hasWildcard = (wildcard in (types1, types2, types3, types4))
                    if match is None or not hasWildcard: # Prefer specific definitions over ones with wildcards
                        match = tordef
                    if not hasWildcard:
                        break
            if match is not None:
                for i in range(len(match.phase)):
                    if match.k[i] != 0:
                        force.addTorsion(torsion[0], torsion[1], torsion[2], torsion[3], match.periodicity[i], match.phase[i], match.k[i])
        for torsion in data.impropers:
            type1 = data.atomType[data.atoms[torsion[0]]]
            type2 = data.atomType[data.atoms[torsion[1]]]
            type3 = data.atomType[data.atoms[torsion[2]]]
            type4 = data.atomType[data.atoms[torsion[3]]]
            done = False
            for tordef in self.improper:
                if done:
                    break
                types1 = tordef.types1
                types2 = tordef.types2
                types3 = tordef.types3
                types4 = tordef.types4
                if type1 in types1:
                    for (t2, t3, t4) in itertools.permutations(((type2, 1), (type3, 2), (type4, 3))):
                        if t2[0] in types2 and t3[0] in types3 and t4[0] in types4:
                            # Workaround to be more consistent with AMBER.  It uses wildcards to define most of its
                            # impropers, which leaves the ordering ambigous.  It then follows some bizarre rules
                            # to pick the order.
                            a1 = torsion[t2[1]]
                            a2 = torsion[t3[1]]
                            e1 = data.atoms[a1].element
                            e2 = data.atoms[a2].element
                            if e1 == e2 and a1 > a2:
                                (a1, a2) = (a2, a1)
                            elif e1 != elem.carbon and (e2 == elem.carbon or e1.mass < e2.mass):
                                (a1, a2) = (a2, a1)
                            for i in range(len(tordef.phase)):
                                if tordef.k[i] != 0:
                                    force.addTorsion(a1, a2, torsion[0], torsion[t4[1]], tordef.periodicity[i], tordef.phase[i], tordef.k[i])
                            done = True
                            break

parsers["PeriodicTorsionForce"] = PeriodicTorsionGenerator.parseElement


class RBTorsion:
    """An RBTorsion records the information for a Ryckaert-Bellemans torsion definition."""

    def __init__(self, types, c):
        self.types1 = types[0]
        self.types2 = types[1]
        self.types3 = types[2]
        self.types4 = types[3]
        self.c = c

class RBTorsionGenerator:
    """An RBTorsionGenerator constructs an RBTorsionForce."""
    
    def __init__(self):
        self.proper = []
        self.improper = []
    
    @staticmethod
    def parseElement(element, ff):
        generator = RBTorsionGenerator()
        generator.ff = ff
        ff._forces.append(generator)
        for torsion in element.findall('Proper'):
            types = ff._findAtomTypes(torsion, 4)
            if types is not None:
                generator.proper.append(RBTorsion(types, [float(torsion.attrib['c'+str(i)]) for i in range(6)]))
        for torsion in element.findall('Improper'):
            types = ff._findAtomTypes(torsion, 4)
            if types is not None:
                generator.improper.append(RBTorsion(types, [float(torsion.attrib['c'+str(i)]) for i in range(6)]))
    
    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
        existing = [sys.getForce(i) for i in range(sys.getNumForces())]
        existing = [f for f in existing if type(f) == mm.RBTorsionForce]
        if len(existing) == 0:
            force = mm.RBTorsionForce()
            sys.addForce(force)
        else:
            force = existing[0]
        wildcard = self.ff._atomClasses['']
        for torsion in data.propers:
            type1 = data.atomType[data.atoms[torsion[0]]]
            type2 = data.atomType[data.atoms[torsion[1]]]
            type3 = data.atomType[data.atoms[torsion[2]]]
            type4 = data.atomType[data.atoms[torsion[3]]]
            match = None
            for tordef in self.proper:
                types1 = tordef.types1
                types2 = tordef.types2
                types3 = tordef.types3
                types4 = tordef.types4
                if (type2 in types2 and type3 in types3 and type4 in types4 and type1 in types1) or (type2 in types3 and type3 in types2 and type4 in types1 and type1 in types4):
                    hasWildcard = (wildcard in (types1, types2, types3, types4))
                    if match is None or not hasWildcard: # Prefer specific definitions over ones with wildcards
                        match = tordef
                    if not hasWildcard:
                        break
            if match is not None:
                force.addTorsion(torsion[0], torsion[1], torsion[2], torsion[3], match.c[0], match.c[1], match.c[2], match.c[3], match.c[4], match.c[5])
        for torsion in data.impropers:
            type1 = data.atomType[data.atoms[torsion[0]]]
            type2 = data.atomType[data.atoms[torsion[1]]]
            type3 = data.atomType[data.atoms[torsion[2]]]
            type4 = data.atomType[data.atoms[torsion[3]]]
            done = False
            for tordef in self.improper:
                if done:
                    break
                types1 = tordef.types1
                types2 = tordef.types2
                types3 = tordef.types3
                types4 = tordef.types4
                if type1 in types1:
                    for (t2, t3, t4) in itertools.permutations(((type2, 1), (type3, 2), (type4, 3))):
                        if t2[0] in types2 and t3[0] in types3 and t4[0] in types4:
                            # Workaround to be more consistent with AMBER.  It uses wildcards to define most of its
                            # impropers, which leaves the ordering ambigous.  It then follows some bizarre rules
                            # to pick the order.
                            a1 = torsion[t2[1]]
                            a2 = torsion[t3[1]]
                            e1 = data.atoms[a1].element
                            e2 = data.atoms[a2].element
                            if e1 == e2 and a1 > a2:
                                (a1, a2) = (a2, a1)
                            elif e1 != elem.carbon and (e2 == elem.carbon or e1.mass < e2.mass):
                                (a1, a2) = (a2, a1)
                            force.addTorsion(a1, a2, torsion[0], torsion[t4[1]], tordef.c[0], tordef.c[1], tordef.c[2], tordef.c[3], tordef.c[4], tordef.c[5])
                            done = True
                            break

parsers["RBTorsionForce"] = RBTorsionGenerator.parseElement


class CMAPTorsion:
    """A CMAPTorsion records the information for a CMAP torsion definition."""

    def __init__(self, types, map):
        self.types1 = types[0]
        self.types2 = types[1]
        self.types3 = types[2]
        self.types4 = types[3]
        self.types5 = types[4]
        self.map = map

class CMAPTorsionGenerator:
    """A CMAPTorsionGenerator constructs a CMAPTorsionForce."""
    
    def __init__(self):
        self.torsions = []
        self.maps = []
    
    @staticmethod
    def parseElement(element, ff):
        generator = CMAPTorsionGenerator()
        generator.ff = ff
        ff._forces.append(generator)
        for map in element.findall('Map'):
            values = [float(x) for x in map.text.split()]
            size = sqrt(len(values))
            if size*size != len(values):
                raise ValueError('CMAP must have the same number of elements along each dimension')
            generator.maps.append(values)
        for torsion in element.findall('Torsion'):
            types = ff._findAtomTypes(torsion, 5)
            if types is not None:
                generator.torsions.append(CMAPTorsion(types, int(torsion.attrib['map'])))
    
    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
        existing = [sys.getForce(i) for i in range(sys.getNumForces())]
        existing = [f for f in existing if type(f) == mm.CMAPTorsionForce]
        if len(existing) == 0:
            force = mm.CMAPTorsionForce()
            sys.addForce(force)
        else:
            force = existing[0]
        for map in self.maps:
            force.addMap(int(sqrt(len(map))), map)
        
        # Find all chains of length 5
        
        uniqueTorsions = set()
        for torsion in data.propers:
            for bond in (data.bonds[x] for x in data.atomBonds[torsion[0]]):
                if bond.atom1 == torsion[0]:
                    atom = bond.atom2
                else:
                    atom = bond.atom1
                if atom != torsion[1]:
                    uniqueTorsions.add((atom, torsion[0], torsion[1], torsion[2], torsion[3]))
            for bond in (data.bonds[x] for x in data.atomBonds[torsion[3]]):
                if bond.atom1 == torsion[3]:
                    atom = bond.atom2
                else:
                    atom = bond.atom1
                if atom != torsion[2]:
                    uniqueTorsions.add((torsion[0], torsion[1], torsion[2], torsion[3], atom))
        torsions = sorted(list(uniqueTorsions))
        wildcard = self.ff._atomClasses['']
        for torsion in torsions:
            type1 = data.atomType[data.atoms[torsion[0]]]
            type2 = data.atomType[data.atoms[torsion[1]]]
            type3 = data.atomType[data.atoms[torsion[2]]]
            type4 = data.atomType[data.atoms[torsion[3]]]
            type5 = data.atomType[data.atoms[torsion[4]]]
            match = None
            for tordef in self.torsions:
                types1 = tordef.types1
                types2 = tordef.types2
                types3 = tordef.types3
                types4 = tordef.types4
                types5 = tordef.types5
                if (type1 in types1 and type2 in types2 and type3 in types3 and type4 in types4 and type5 in types5) or (type1 in types5 and type2 in types4 and type3 in types3 and type4 in types2 and type5 in types1):
                    hasWildcard = (wildcard in (types1, types2, types3, types4, types5))
                    if match is None or not hasWildcard: # Prefer specific definitions over ones with wildcards
                        match = tordef
                    if not hasWildcard:
                        break
            if match is not None:
                force.addTorsion(match.map, torsion[0], torsion[1], torsion[2], torsion[3], torsion[1], torsion[2], torsion[3], torsion[4])

parsers["CMAPTorsionForce"] = CMAPTorsionGenerator.parseElement


class NonbondedGenerator:
    """A NonbondedGenerator constructs a NonbondedForce."""
    
    def __init__(self, coulomb14scale, lj14scale):
        self.coulomb14scale = coulomb14scale
        self.lj14scale = lj14scale
        self.typeMap = {}

    @staticmethod
    def parseElement(element, ff):
        existing = [f for f in ff._forces if isinstance(f, NonbondedGenerator)]
        if len(existing) == 0:
            generator = NonbondedGenerator(float(element.attrib['coulomb14scale']), float(element.attrib['lj14scale']))
            ff._forces.append(generator)
        else:
            # Multiple <NonbondedForce> tags were found, probably in different files.  Simply add more types to the existing one.
            generator = existing[0]
            if generator.coulomb14scale != float(element.attrib['coulomb14scale']) or generator.lj14scale != float(element.attrib['lj14scale']):
                raise ValueError('Found multiple NonbondedForce tags with different 1-4 scales') 
        for atom in element.findall('Atom'):
            types = ff._findAtomTypes(atom, 1)
            if types is not None:
                values = (float(atom.attrib['charge']), float(atom.attrib['sigma']), float(atom.attrib['epsilon']))
                for t in types[0]:
                    generator.typeMap[t] = values
    
    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
        methodMap = {NoCutoff:mm.NonbondedForce.NoCutoff,
                     CutoffNonPeriodic:mm.NonbondedForce.CutoffNonPeriodic,
                     CutoffPeriodic:mm.NonbondedForce.CutoffPeriodic,
                     Ewald:mm.NonbondedForce.Ewald,
                     PME:mm.NonbondedForce.PME}
        if nonbondedMethod not in methodMap:
            raise ValueError('Illegal nonbonded method for NonbondedForce') 
        force = mm.NonbondedForce()
        for atom in data.atoms:
            t = data.atomType[atom]
            if t in self.typeMap:
                values = self.typeMap[t]
                force.addParticle(values[0], values[1], values[2])
            else:
                raise ValueError('No nonbonded parameters defined for atom type '+t) 
        bondIndices = []
        for bond in data.bonds:
            bondIndices.append((bond.atom1, bond.atom2))
        force.createExceptionsFromBonds(bondIndices, self.coulomb14scale, self.lj14scale)
        force.setNonbondedMethod(methodMap[nonbondedMethod])
        force.setCutoffDistance(nonbondedCutoff)
        if 'ewaldErrorTolerance' in args:
            force.setEwaldErrorTolerance(args['ewaldErrorTolerance'])
        sys.addForce(force)

parsers["NonbondedForce"] = NonbondedGenerator.parseElement


class GBSAOBCGenerator:
    """A GBSAOBCGenerator constructs a GBSAOBCForce."""
    
    def __init__(self):
        self.typeMap = {}

    @staticmethod
    def parseElement(element, ff):
        generator = GBSAOBCGenerator()
        ff._forces.append(generator)
        for atom in element.findall('Atom'):
            types = ff._findAtomTypes(atom, 1)
            if types is not None:
                values = (float(atom.attrib['charge']), float(atom.attrib['radius']), float(atom.attrib['scale']))
                for t in types[0]:
                    generator.typeMap[t] = values
    
    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
        methodMap = {NoCutoff:mm.NonbondedForce.NoCutoff,
                     CutoffNonPeriodic:mm.NonbondedForce.CutoffNonPeriodic,
                     CutoffPeriodic:mm.NonbondedForce.CutoffPeriodic}
        if nonbondedMethod not in methodMap:
            raise ValueError('Illegal nonbonded method for GBSAOBCForce') 
        force = mm.GBSAOBCForce()
        for atom in data.atoms:
            t = data.atomType[atom]
            if t in self.typeMap:
                values = self.typeMap[t]
                force.addParticle(values[0], values[1], values[2])
            else:
                raise ValueError('No GBSAOBC parameters defined for atom type '+t) 
        force.setNonbondedMethod(methodMap[nonbondedMethod])
        force.setCutoffDistance(nonbondedCutoff)
        sys.addForce(force)

parsers["GBSAOBCForce"] = GBSAOBCGenerator.parseElement


class GBVIGenerator:

    """A GBVIGenerator constructs a GBVIForce."""
    
    def __init__(self,ff):

Peter Eastman's avatar
Peter Eastman committed
926
927
928
929
930
        self.ff = ff
        self.fixedParameters = {}
        self.fixedParameters['soluteDielectric'] = 1.0
        self.fixedParameters['solventDielectric'] = 78.3
        self.fixedParameters['scalingMethod'] = 1
931
        self.fixedParameters['quinticUpperBornRadiusLimit'] = 5.0
Peter Eastman's avatar
Peter Eastman committed
932
        self.fixedParameters['quinticLowerLimitFactor'] = 0.8
933

Peter Eastman's avatar
Peter Eastman committed
934
        self.typeMap = {}
935
936
937

    @staticmethod
    def parseElement(element, ff):
Peter Eastman's avatar
Peter Eastman committed
938
        generator = GBVIGenerator(ff)
939
        for key in generator.fixedParameters.iterkeys():
Peter Eastman's avatar
Peter Eastman committed
940
941
            if (key in element.attrib):
                generator.fixedParameters[key] = float(element.attrib[key])
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974

        ff._forces.append(generator)
        for atom in element.findall('Atom'):
            types = ff._findAtomTypes(atom, 1)
            if types is not None:
                values = (float(atom.attrib['charge']), float(atom.attrib['radius']), float(atom.attrib['gamma']))
                for t in types[0]:
                    generator.typeMap[t] = values
    
    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):

        methodMap = {NoCutoff:mm.NonbondedForce.NoCutoff,
                     CutoffNonPeriodic:mm.NonbondedForce.CutoffNonPeriodic,
                     CutoffPeriodic:mm.NonbondedForce.CutoffPeriodic}

        if nonbondedMethod not in methodMap:
            raise ValueError('Illegal nonbonded method for GB/VI Force') 

        # add particles

        force = mm.GBVIForce()
        for atom in data.atoms:
            t = data.atomType[atom]
            if t in self.typeMap:
                values = self.typeMap[t]
                force.addParticle(values[0], values[1], values[2])
            else:
                raise ValueError('No GB/VI parameters defined for atom type '+t) 

        # get HarmonicBond generator -- exit if not found

        hbGenerator = 0
        for generator in self.ff._forces:
Peter Eastman's avatar
Peter Eastman committed
975
            if (generator.__class__.__name__ == 'HarmonicBondGenerator'): 
976
977
978
               hbGenerator = generator
               break

Peter Eastman's avatar
Peter Eastman committed
979
        if (hbGenerator == 0):
980
981
982
983
984
985
986
            raise ValueError('HarmonicBondGenerator not found.') 

        # add bonds

        for bond in data.bonds:
            type1 = data.atomType[data.atoms[bond.atom1]]
            type2 = data.atomType[data.atoms[bond.atom2]]
Peter Eastman's avatar
Peter Eastman committed
987
            hit = 0
988
989
990
991
992
993
994
995
996
            for i in range(len(hbGenerator.types1)):
                types1 = hbGenerator.types1[i]
                types2 = hbGenerator.types2[i]
                if (type1 in types1 and type2 in types2) or (type1 in types2 and type2 in types1):
                    #bond.length = hbGenerator.length[i]
                    force.addBond(bond.atom1, bond.atom2, hbGenerator.length[i])

        force.setNonbondedMethod(methodMap[nonbondedMethod])
        force.setCutoffDistance(nonbondedCutoff)
Peter Eastman's avatar
Peter Eastman committed
997
998
999
1000
1001
        force.setSolventDielectric(self.fixedParameters['solventDielectric'])
        force.setSoluteDielectric(self.fixedParameters['soluteDielectric'])
        force.setBornRadiusScalingMethod(self.fixedParameters['scalingMethod'])
        force.setQuinticLowerLimitFactor(self.fixedParameters['quinticLowerLimitFactor'])
        force.setQuinticUpperBornRadiusLimit(self.fixedParameters['quinticUpperBornRadiusLimit'])
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
        
        sys.addForce(force)

parsers["GBVIForce"] = GBVIGenerator.parseElement

class CustomBondGenerator:
    """A CustomBondGenerator constructs a CustomBondForce."""
    
    def __init__(self):
        self.types1 = []
        self.types2 = []
        self.globalParams = {}
        self.perBondParams = []
        self.paramValues = []
    
    @staticmethod
    def parseElement(element, ff):
        generator = CustomBondGenerator()
        ff._forces.append(generator)
        generator.energy = element.attrib['energy']
        for param in element.findall('GlobalParameter'):
            generator.globalParams[param.attrib['name']] = float(param.attrib['defaultValue'])
        for param in element.findall('PerBondParameter'):
            generator.perBondParams.append(param.attrib['name'])
        for bond in element.findall('Bond'):
            types = ff._findAtomTypes(bond, 2)
            if types is not None:
                generator.types1.append(types[0])
                generator.types2.append(types[1])
                generator.paramValues.append([float(bond.attrib[param]) for param in generator.perBondParams])
    
    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
        force = mm.CustomBondForce(self.energy)
        sys.addForce(force)
        for param in self.globalParams:
            force.addGlobalParameter(param, self.globalParams[param])
        for param in self.perBondParams:
            force.addPerBondParameter(param)
        for bond in data.bonds:
            type1 = data.atomType[data.atoms[bond.atom1]]
            type2 = data.atomType[data.atoms[bond.atom2]]
            for i in range(len(self.types1)):
                types1 = self.types1[i]
                types2 = self.types2[i]
                if (type1 in types1 and type2 in types2) or (type1 in types2 and type2 in types1):
                    force.addBond(bond.atom1, bond.atom2, self.paramValues[i])
                    break

parsers["CustomBondForce"] = CustomBondGenerator.parseElement


class CustomAngleGenerator:
    """A CustomAngleGenerator constructs a CustomAngleForce."""
    
    def __init__(self):
        self.types1 = []
        self.types2 = []
        self.types3 = []
        self.globalParams = {}
        self.perAngleParams = []
        self.paramValues = []
    
    @staticmethod
    def parseElement(element, ff):
        generator = CustomAngleGenerator()
        ff._forces.append(generator)
        generator.energy = element.attrib['energy']
        for param in element.findall('GlobalParameter'):
            generator.globalParams[param.attrib['name']] = float(param.attrib['defaultValue'])
        for param in element.findall('PerAngleParameter'):
            generator.perAngleParams.append(param.attrib['name'])
        for angle in element.findall('Angle'):
            types = ff._findAtomTypes(angle, 3)
            if types is not None:
                generator.types1.append(types[0])
                generator.types2.append(types[1])
                generator.types3.append(types[2])
                generator.paramValues.append([float(angle.attrib[param]) for param in generator.perAngleParams])
    
    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
        force = mm.CustomAngleForce(self.energy)
        sys.addForce(force)
        for param in self.globalParams:
            force.addGlobalParameter(param, self.globalParams[param])
        for param in self.perAngleParams:
            force.addPerAngleParameter(param)
        for angle in data.angles:
            type1 = data.atomType[data.atoms[angle[0]]]
            type2 = data.atomType[data.atoms[angle[1]]]
            type3 = data.atomType[data.atoms[angle[2]]]
            for i in range(len(self.types1)):
                types1 = self.types1[i]
                types2 = self.types2[i]
                types3 = self.types3[i]
                if (type1 in types1 and type2 in types2 and type3 in types3) or (type1 in types3 and type2 in types2 and type3 in types1):
                    force.addAngle(angle[0], angle[1], angle[2], self.paramValues[i])
                    break

parsers["CustomAngleForce"] = CustomAngleGenerator.parseElement


class CustomTorsion:
    """A CustomTorsion records the information for a custom torsion definition."""

    def __init__(self, types, paramValues):
        self.types1 = types[0]
        self.types2 = types[1]
        self.types3 = types[2]
        self.types4 = types[3]
        self.paramValues = paramValues

class CustomTorsionGenerator:
    """A CustomTorsionGenerator constructs a CustomTorsionForce."""
    
    def __init__(self):
        self.proper = []
        self.improper = []
        self.globalParams = {}
        self.perTorsionParams = []
    
    @staticmethod
    def parseElement(element, ff):
        generator = CustomTorsionGenerator()
        generator.ff = ff
        ff._forces.append(generator)
        generator.energy = element.attrib['energy']
        for param in element.findall('GlobalParameter'):
            generator.globalParams[param.attrib['name']] = float(param.attrib['defaultValue'])
        for param in element.findall('PerTorsionParameter'):
            generator.perTorsionParams.append(param.attrib['name'])
        for torsion in element.findall('Proper'):
            types = ff._findAtomTypes(torsion, 4)
            if types is not None:
                generator.proper.append(CustomTorsion(types, [float(torsion.attrib[param]) for param in generator.perTorsionParams]))
        for torsion in element.findall('Improper'):
            types = ff._findAtomTypes(torsion, 4)
            if types is not None:
                generator.improper.append(CustomTorsion(types, [float(torsion.attrib[param]) for param in generator.perTorsionParams]))
    
    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
        force = mm.CustomTorsionForce(self.energy)
        sys.addForce(force)
        for param in self.globalParams:
            force.addGlobalParameter(param, self.globalParams[param])
        for param in self.perTorsionParams:
            force.addPerTorsionParameter(param)
        wildcard = self.ff._atomClasses['']
        for torsion in data.propers:
            type1 = data.atomType[data.atoms[torsion[0]]]
            type2 = data.atomType[data.atoms[torsion[1]]]
            type3 = data.atomType[data.atoms[torsion[2]]]
            type4 = data.atomType[data.atoms[torsion[3]]]
            match = None
            for tordef in self.proper:
                types1 = tordef.types1
                types2 = tordef.types2
                types3 = tordef.types3
                types4 = tordef.types4
                if (type2 in types2 and type3 in types3 and type4 in types4 and type1 in types1) or (type2 in types3 and type3 in types2 and type4 in types1 and type1 in types4):
                    hasWildcard = (wildcard in (types1, types2, types3, types4))
                    if match is None or not hasWildcard: # Prefer specific definitions over ones with wildcards
                        match = tordef
                    if not hasWildcard:
                        break
            if match is not None:
                force.addTorsion(torsion[0], torsion[1], torsion[2], torsion[3], match.paramValues)
        for torsion in data.impropers:
            type1 = data.atomType[data.atoms[torsion[0]]]
            type2 = data.atomType[data.atoms[torsion[1]]]
            type3 = data.atomType[data.atoms[torsion[2]]]
            type4 = data.atomType[data.atoms[torsion[3]]]
            done = False
            for tordef in self.improper:
                if done:
                    break
                types1 = tordef.types1
                types2 = tordef.types2
                types3 = tordef.types3
                types4 = tordef.types4
                if type1 in types1:
                    for (t2, t3, t4) in itertools.permutations(((type2, 1), (type3, 2), (type4, 3))):
                        if t2[0] in types2 and t3[0] in types3 and t4[0] in types4:
                            # Workaround to be more consistent with AMBER.  It uses wildcards to define most of its
                            # impropers, which leaves the ordering ambigous.  It then follows some bizarre rules
                            # to pick the order.
                            a1 = torsion[t2[1]]
                            a2 = torsion[t3[1]]
                            e1 = data.atoms[a1].element
                            e2 = data.atoms[a2].element
                            if e1 == e2 and a1 > a2:
                                (a1, a2) = (a2, a1)
                            elif e1 != elem.carbon and (e2 == elem.carbon or e1.mass < e2.mass):
                                (a1, a2) = (a2, a1)
                            force.addTorsion(a1, a2, torsion[0], torsion[t4[1]], tordef.paramValues)
                            done = True
                            break

parsers["CustomTorsionForce"] = CustomTorsionGenerator.parseElement


class CustomGBGenerator:
    """A CustomGBGenerator constructs a CustomGBForce."""
    
    def __init__(self):
        self.typeMap = {}
        self.globalParams = {}
        self.perParticleParams = []
        self.paramValues = []
        self.computedValues = []
        self.energyTerms = []
        self.functions = []

    @staticmethod
    def parseElement(element, ff):
        generator = CustomGBGenerator()
        ff._forces.append(generator)
        for param in element.findall('GlobalParameter'):
            generator.globalParams[param.attrib['name']] = float(param.attrib['defaultValue'])
        for param in element.findall('PerParticleParameter'):
            generator.perParticleParams.append(param.attrib['name'])
        for atom in element.findall('Atom'):
            types = ff._findAtomTypes(atom, 1)
            if types is not None:
                values = [float(atom.attrib[param]) for param in generator.perParticleParams]
                for t in types[0]:
                    generator.typeMap[t] = values
        computationMap = {"SingleParticle" : mm.CustomGBForce.SingleParticle,
                          "ParticlePair" : mm.CustomGBForce.ParticlePair,
                          "ParticlePairNoExclusions" : mm.CustomGBForce.ParticlePairNoExclusions}
        for value in element.findall('ComputedValue'):
            generator.computedValues.append((value.attrib['name'], value.text, computationMap[value.attrib['type']]))
        for term in element.findall('EnergyTerm'):
            generator.energyTerms.append((term.text, computationMap[term.attrib['type']]))
        for function in element.findall("Function"):
            values = [float(x) for x in function.text.split()]
            generator.functions.append((function.attrib['name'], values, float(function.attrib['min']), float(function.attrib['max'])))
    
    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
        methodMap = {NoCutoff:mm.CustomGBForce.NoCutoff,
                     CutoffNonPeriodic:mm.CustomGBForce.CutoffNonPeriodic,
                     CutoffPeriodic:mm.CustomGBForce.CutoffPeriodic}
        if nonbondedMethod not in methodMap:
            raise ValueError('Illegal nonbonded method for CustomGBForce') 
        force = mm.CustomGBForce()
        for param in self.globalParams:
            force.addGlobalParameter(param, self.globalParams[param])
        for param in self.perParticleParams:
            force.addPerParticleParameter(param)
        for value in self.computedValues:
            force.addComputedValue(value[0], value[1], value[2])
        for term in self.energyTerms:
            force.addEnergyTerm(term[0], term[1])
        for function in self.functions:
            force.addFunction(function[0], function[1], function[2], function[3])
        for atom in data.atoms:
            t = data.atomType[atom]
            if t in self.typeMap:
                values = self.typeMap[t]
                force.addParticle(self.typeMap[t])
            else:
                raise ValueError('No CustomGB parameters defined for atom type '+t) 
        force.setNonbondedMethod(methodMap[nonbondedMethod])
        force.setCutoffDistance(nonbondedCutoff)
        sys.addForce(force)

parsers["CustomGBForce"] = CustomGBGenerator.parseElement

Peter Eastman's avatar
Peter Eastman committed
1269
def getAtomPrint(data, atomIndex):
1270

Peter Eastman's avatar
Peter Eastman committed
1271
1272
1273
    if (atomIndex < len(data.atoms)):
        atom = data.atoms[atomIndex]
        returnString = "%4s %4s %5d" % (atom.name, atom.residue.name, atom.residue.index)
1274
    else:
Peter Eastman's avatar
Peter Eastman committed
1275
        returnString = "NA"
1276
1277
1278
1279
1280

    return returnString

#=============================================================================================

Peter Eastman's avatar
Peter Eastman committed
1281
def countConstraint(data):
1282

Peter Eastman's avatar
Peter Eastman committed
1283
    bondCount = 0
1284
1285
1286
1287
1288
1289
1290
    angleCount = 0
    for bond in data.bonds:
        if bond.isConstrained:
            bondCount += 1

    angleCount = 0
    for (angle, isConstrained) in zip(data.angles, data.isAngleConstrained):
Peter Eastman's avatar
Peter Eastman committed
1291
        if (isConstrained):
1292
1293
            angleCount += 1
 
Peter Eastman's avatar
Peter Eastman committed
1294
    print "Constraints bond=%d angle=%d  total=%d" % (bondCount, angleCount, (bondCount+angleCount))
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305

class AmoebaHarmonicBondGenerator:

    #=============================================================================================

    """An AmoebaHarmonicBondGenerator constructs a AmoebaHarmonicBondForce."""

    #=============================================================================================
    
    def __init__(self, cubic, quartic):

Peter Eastman's avatar
Peter Eastman committed
1306
1307
1308
1309
1310
1311
        self.cubic = cubic
        self.quartic = quartic
        self.types1 = []
        self.types2 = []
        self.length = []
        self.k = []
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
        self.hasBeenCalled = 0
    
    #=============================================================================================

    @staticmethod
    def parseElement(element, forceField):

        # <AmoebaHarmonicBondForce bond-cubic="-25.5" bond-quartic="379.3125">
        # <Bond class1="1" class2="2" length="0.1437" k="156900.0"/>
    
Peter Eastman's avatar
Peter Eastman committed
1322
        generator = AmoebaHarmonicBondGenerator(float(element.attrib['bond-cubic']), float(element.attrib['bond-quartic']))
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
        forceField._forces.append(generator)
        for bond in element.findall('Bond'):
            types = forceField._findAtomTypes(bond, 2)
            if types is not None:
                generator.types1.append(types[0])
                generator.types2.append(types[1])
                generator.length.append(float(bond.attrib['length']))
                generator.k.append(float(bond.attrib['k']))
            else:
                outputString = "AmoebaHarmonicBondGenerator: error getting types: %s %s" % (
                                    bond.attrib['class1'],
Peter Eastman's avatar
Peter Eastman committed
1334
1335
                                    bond.attrib['class2'])
                raise ValueError(outputString) 
1336
1337
1338
    
    #=============================================================================================

Peter Eastman's avatar
Peter Eastman committed
1339
    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
1340

Peter Eastman's avatar
Peter Eastman committed
1341
1342
        verbose = 0
        if (self.hasBeenCalled):
1343
1344
             return

Peter Eastman's avatar
Peter Eastman committed
1345
1346
        if (verbose): 
            countConstraint(data)
1347
1348
1349

        self.hasBeenCalled = 1

Peter Eastman's avatar
Peter Eastman committed
1350
1351
        existing = [sys.getForce(i) for i in range(sys.getNumForces())]
        existing = [f for f in existing if type(f) == mm.AmoebaHarmonicBondForce]
1352
1353
1354
1355
1356
1357
        if len(existing) == 0:
            force = mm.AmoebaHarmonicBondForce()
            sys.addForce(force)
        else:
            force = existing[0]

Peter Eastman's avatar
Peter Eastman committed
1358
1359
        force.setAmoebaGlobalHarmonicBondCubic(self.cubic)
        force.setAmoebaGlobalHarmonicBondQuartic(self.quartic)
1360

Peter Eastman's avatar
Peter Eastman committed
1361
        if (verbose):
1362
1363
1364
1365
1366
1367
            print "In AmoebaHarmonicBondGenerator bonds=%d " % (len(data.bonds))

        count = 0
        for bond in data.bonds:
            type1 = data.atomType[data.atoms[bond.atom1]]
            type2 = data.atomType[data.atoms[bond.atom2]]
Peter Eastman's avatar
Peter Eastman committed
1368
            hit = 0
1369
1370
1371
1372
1373
            for i in range(len(self.types1)):
                types1 = self.types1[i]
                types2 = self.types2[i]
                if (type1 in types1 and type2 in types2) or (type1 in types2 and type2 in types1):
                    bond.length = self.length[i]
Peter Eastman's avatar
Peter Eastman committed
1374
                    hit = 1
1375
1376
                    if bond.isConstrained:
                        sys.addConstraint(bond.atom1, bond.atom2, self.length[i])
Peter Eastman's avatar
Peter Eastman committed
1377
1378
1379
                        if (verbose):
                            atomS1 = getAtomPrint(data, bond.atom1)
                            atomS2 = getAtomPrint(data, bond.atom2)
1380
1381
                            print "AmoebaHarmonicBondGenerator %5d %5d %5d [%s %s] [%5s %5s] %15.6f %15.6f Constraint" % (count, bond.atom1, bond.atom2, atomS1, atomS2, type1, type2, self.length[i], self.k[i])
                    elif self.k[i] != 0:
Peter Eastman's avatar
Peter Eastman committed
1382
1383
1384
                        if (verbose):
                            atomS1 = getAtomPrint(data, bond.atom1)
                            atomS2 = getAtomPrint(data, bond.atom2)
1385
1386
1387
                            print "AmoebaHarmonicBondGenerator %5d %5d %5d [%s %s] [%5s %5s] %15.6f %15.6f" % (count, bond.atom1, bond.atom2, atomS1, atomS2, type1, type2, self.length[i], self.k[i])
                        force.addBond(bond.atom1, bond.atom2, self.length[i], self.k[i])
                    break
Peter Eastman's avatar
Peter Eastman committed
1388
            if (hit == 0): 
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
                print "AmoebaHarmonicBondGenerator missing: %5d types=[%5s %5s] atoms=[%6d %6d] " % (count, type1, type2, bond.atom1, bond.atom2)

            count += 1

parsers["AmoebaHarmonicBondForce"] = AmoebaHarmonicBondGenerator.parseElement

#=============================================================================================
# Add angle constraint
#=============================================================================================
    
Peter Eastman's avatar
Peter Eastman committed
1399
def addAngleConstraint(angle, idealAngle, data, sys):
1400
1401
1402

    # Find the two bonds that make this angle.
                    
Peter Eastman's avatar
Peter Eastman committed
1403
1404
    bond1 = None
    bond2 = None
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
    for bond in data.atomBonds[angle[1]]:
        atom1 = data.bonds[bond].atom1
        atom2 = data.bonds[bond].atom2
        if atom1 == angle[0] or atom2 == angle[0]:
            bond1 = bond
        elif atom1 == angle[2] or atom2 == angle[2]:
            bond2 = bond
                    
        # Compute the distance between atoms and add a constraint
                    
        if bond1 is not None and bond2 is not None:
            l1 = data.bonds[bond1].length
            l2 = data.bonds[bond2].length
            if l1 is not None and l2 is not None:
Peter Eastman's avatar
Peter Eastman committed
1419
                length = sqrt(l1*l1 + l2*l2 - 2*l1*l2*cos(idealAngle))
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432
                sys.addConstraint(angle[0], angle[2], length)
                return


#=============================================================================================
class AmoebaHarmonicAngleGenerator:

    #=============================================================================================
    """An AmoebaHarmonicAngleGenerator constructs a AmoebaHarmonicAngleForce."""
    #=============================================================================================
    
    def __init__(self, forceField, cubic, quartic, pentic, sextic):

Peter Eastman's avatar
Peter Eastman committed
1433
1434
1435
1436
1437
        self.forceField = forceField
        self.cubic = cubic
        self.quartic = quartic
        self.pentic = pentic
        self.sextic = sextic
1438

Peter Eastman's avatar
Peter Eastman committed
1439
1440
1441
        self.types1 = []
        self.types2 = []
        self.types3 = []
1442

Peter Eastman's avatar
Peter Eastman committed
1443
1444
        self.angle = []
        self.k = []
1445
1446
1447
1448
1449
1450
1451
1452
1453
1454
1455

        self.hasBeenCalled = 0
    
    #=============================================================================================

    @staticmethod
    def parseElement(element, forceField):

        # <AmoebaHarmonicAngleForce angle-cubic="-0.014" angle-quartic="5.6e-05" angle-pentic="-7e-07" angle-sextic="2.2e-08">
        #   <Angle class1="2" class2="1" class3="3" k="0.0637259642196" angle1="122.00"  />

Peter Eastman's avatar
Peter Eastman committed
1456
        generator = AmoebaHarmonicAngleGenerator(forceField, float(element.attrib['angle-cubic']), float(element.attrib['angle-quartic']),  float(element.attrib['angle-pentic']), float(element.attrib['angle-sextic']))
1457
1458
1459
1460
1461
1462
1463
1464
1465
1466
        forceField._forces.append(generator)
        for angle in element.findall('Angle'):
            types = forceField._findAtomTypes(angle, 3)
            if types is not None:

                generator.types1.append(types[0])
                generator.types2.append(types[1])
                generator.types3.append(types[2])

                angleList = []
Peter Eastman's avatar
Peter Eastman committed
1467
                angleList.append(float(angle.attrib['angle1']))
1468
1469

                try:
Peter Eastman's avatar
Peter Eastman committed
1470
                    angleList.append(float(angle.attrib['angle2']))
1471
                    try:
Peter Eastman's avatar
Peter Eastman committed
1472
                        angleList.append(float(angle.attrib['angle3']))
1473
1474
1475
1476
1477
1478
1479
1480
1481
1482
                    except:
                        pass
                except:
                    pass
                generator.angle.append(angleList)
                generator.k.append(float(angle.attrib['k']))
            else:
                outputString = "AmoebaHarmonicAngleGenerator: error getting types: %s %s %s" % (
                                    angle.attrib['class1'],
                                    angle.attrib['class2'],
Peter Eastman's avatar
Peter Eastman committed
1483
1484
                                    angle.attrib['class3'])
                raise ValueError(outputString) 
1485
1486
1487
1488
1489
1490
    
    #=============================================================================================
    # createForce is bypassed here since the AmoebaOutOfPlaneBendForce generator must first execute
    # and partition angles into in-plane and non-in-plane angles
    #=============================================================================================
    
Peter Eastman's avatar
Peter Eastman committed
1491
    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
1492
1493
1494
1495
1496
1497
1498
        pass

    #=============================================================================================
    # createForcePostOpBendAngle is called by AmoebaOutOfPlaneBendForce with the list of
    # non-in-plane angles
    #=============================================================================================
    
Peter Eastman's avatar
Peter Eastman committed
1499
    def createForcePostOpBendAngle(self, sys, data, nonbondedMethod, nonbondedCutoff, angleList, args):
1500

Peter Eastman's avatar
Peter Eastman committed
1501
1502
        verbose = 0
        if (self.hasBeenCalled):
1503
1504
1505
1506
1507
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
1518
             return
        self.hasBeenCalled += 1

        # get force

        existing = [sys.getForce(i) for i in range(sys.getNumForces())]
        existing = [f for f in existing if type(f) == mm.AmoebaHarmonicAngleForce]

        if len(existing) == 0:
            force = mm.AmoebaHarmonicAngleForce()
            sys.addForce(force)
        else:
            force = existing[0]

        # scalars

Peter Eastman's avatar
Peter Eastman committed
1519
1520
1521
1522
        force.setAmoebaGlobalHarmonicAngleCubic(self.cubic)
        force.setAmoebaGlobalHarmonicAngleQuartic(self.quartic)
        force.setAmoebaGlobalHarmonicAnglePentic(self.pentic)
        force.setAmoebaGlobalHarmonicAngleSextic(self.sextic)
1523

Peter Eastman's avatar
Peter Eastman committed
1524
        if (verbose):
1525
1526
1527
1528
            print "In AmoebaHarmonicAngleGenerator angles=%d " % (len(data.angles))

        count = 0
        for angleDict in angleList:
Peter Eastman's avatar
Peter Eastman committed
1529
1530
            angle = angleDict['angle']
            isConstrained = angleDict['isConstrained']
1531

Peter Eastman's avatar
Peter Eastman committed
1532
1533
1534
1535
            type1 = data.atomType[data.atoms[angle[0]]]
            type2 = data.atomType[data.atoms[angle[1]]]
            type3 = data.atomType[data.atoms[angle[2]]]
            hit = 0
1536
1537
1538
1539
1540
            for i in range(len(self.types1)):
                types1 = self.types1[i]
                types2 = self.types2[i]
                types3 = self.types3[i]
                if (type1 in types1 and type2 in types2 and type3 in types3) or (type1 in types3 and type2 in types2 and type3 in types1):
Peter Eastman's avatar
Peter Eastman committed
1541
                    hit = 1
1542
1543
                    if isConstrained and self.k[i] != 0.0:
                        angleDict['idealAngle'] = self.angle[i][0]
Peter Eastman's avatar
Peter Eastman committed
1544
                        addAngleConstraint(angle, self.angle[i][0], data, sys)
1545
                    elif self.k[i] != 0:
Peter Eastman's avatar
Peter Eastman committed
1546
1547
                        lenAngle = len(self.angle[i])
                        if (lenAngle > 1):
1548
1549
1550
1551
1552
1553
                            # get k-index by counting number of non-angle hydrogens on the central atom
                            # based on kangle.f
                            numberOfHydrogens = 0
                            for bond in data.atomBonds[angle[1]]:
                                atom1 = data.bonds[bond].atom1
                                atom2 = data.bonds[bond].atom2
Peter Eastman's avatar
Peter Eastman committed
1554
                                if (atom1 == angle[1] and atom2 != angle[0] and atom2 != angle[2] and (sys.getParticleMass(atom2)/unit.dalton) < 1.90):
1555
                                    numberOfHydrogens += 1
Peter Eastman's avatar
Peter Eastman committed
1556
                                if (atom2 == angle[1] and atom1 != angle[0] and atom1 != angle[2] and (sys.getParticleMass(atom1)/unit.dalton) < 1.90):
1557
                                    numberOfHydrogens += 1
Peter Eastman's avatar
Peter Eastman committed
1558
                            if (numberOfHydrogens < lenAngle):
1559
1560
1561
1562
1563
1564
1565
1566
                                angleValue =  self.angle[i][numberOfHydrogens]
                            else:
                                print "Error: AmoebaHarmonicAngleGenerator angle index=%d is out of range: [0, %5d] " % (numberOfHydrogens, lenAngle)
                                sys.exit(-1)
                        else:
                            angleValue =  self.angle[i][0]
               
                        angleDict['idealAngle'] = angleValue
Peter Eastman's avatar
Peter Eastman committed
1567
                        force.addAngle(angle[0], angle[1], angle[2], angleValue, self.k[i])
1568
                    break
Peter Eastman's avatar
Peter Eastman committed
1569
            if (hit == 0): 
1570
1571
1572
1573
1574
1575
1576
1577
1578
                print "AmoebaHarmonicAngleGenerator missing: %5d %s %s %s " % (count, type1, type2, type3)

            count += 1

    #=============================================================================================
    # createForcePostOpBendInPlaneAngle is called by AmoebaOutOfPlaneBendForce with the list of
    # in-plane angles
    #=============================================================================================
    
Peter Eastman's avatar
Peter Eastman committed
1579
    def createForcePostOpBendInPlaneAngle(self, sys, data, nonbondedMethod, nonbondedCutoff, angleList, args):
1580

Peter Eastman's avatar
Peter Eastman committed
1581
        verbose = 0
1582
1583
1584
1585
1586
1587
1588
1589
1590
1591
1592
1593
1594
1595
1596
        self.hasBeenCalled += 1

        # get force

        existing = [sys.getForce(i) for i in range(sys.getNumForces())]
        existing = [f for f in existing if type(f) == mm.AmoebaHarmonicInPlaneAngleForce]

        if len(existing) == 0:
            force = mm.AmoebaHarmonicInPlaneAngleForce()
            sys.addForce(force)
        else:
            force = existing[0]

        # scalars

Peter Eastman's avatar
Peter Eastman committed
1597
1598
1599
1600
        force.setAmoebaGlobalHarmonicInPlaneAngleCubic(self.cubic)
        force.setAmoebaGlobalHarmonicInPlaneAngleQuartic(self.quartic)
        force.setAmoebaGlobalHarmonicInPlaneAnglePentic(self.pentic)
        force.setAmoebaGlobalHarmonicInPlaneAngleSextic(self.sextic)
1601

Peter Eastman's avatar
Peter Eastman committed
1602
        if (verbose):
1603
1604
1605
1606
1607
            print "In AmoebaHarmonicAngleGenerator angles=%d " % (len(data.angles))

        count = 0
        for angleDict in angleList:
 
Peter Eastman's avatar
Peter Eastman committed
1608
1609
            angle = angleDict['angle']
            isConstrained = angleDict['isConstrained']
1610

Peter Eastman's avatar
Peter Eastman committed
1611
1612
1613
            type1 = data.atomType[data.atoms[angle[0]]]
            type2 = data.atomType[data.atoms[angle[1]]]
            type3 = data.atomType[data.atoms[angle[2]]]
1614

Peter Eastman's avatar
Peter Eastman committed
1615
            hit = 0
1616
1617
1618
1619
1620
1621
1622
            for i in range(len(self.types1)):

                types1 = self.types1[i]
                types2 = self.types2[i]
                types3 = self.types3[i]

                if (type1 in types1 and type2 in types2 and type3 in types3) or (type1 in types3 and type2 in types2 and type3 in types1):
Peter Eastman's avatar
Peter Eastman committed
1623
                    if (verbose):
1624
                        print "AmoebaHarmonicInPlaneAngleGenerator %5d %5d %5d %5d len=%d [%5s %5s %5s] %15.6f %15.6f" % (count, angle[0], angle[1], angle[2], len(angle), type1, type2, type3, self.angle[i][0], self.k[i])
Peter Eastman's avatar
Peter Eastman committed
1625
                    hit = 1
1626
                    angleDict['idealAngle'] = self.angle[i][0]
Peter Eastman's avatar
Peter Eastman committed
1627
1628
                    if (isConstrained and self.k[i] != 0.0):
                        addAngleConstraint(angle, self.angle[i][0], data, sys)
1629
1630
1631
1632
                    else:
                        force.addAngle(angle[0], angle[1], angle[2], angle[3], self.angle[i][0], self.k[i])
                    break

Peter Eastman's avatar
Peter Eastman committed
1633
            if (hit == 0): 
1634
1635
1636
1637
1638
1639
1640
1641
1642
1643
1644
1645
1646
1647
1648
1649
1650
1651
1652
1653
1654
1655
                print "AmoebaHarmonicInPlaneAngleGenerator missing: %5d %s %s " % (count, type1, type2, type3)

            count += 1

parsers["AmoebaHarmonicAngleForce"] = AmoebaHarmonicAngleGenerator.parseElement

#=============================================================================================
# Generator for the AmoebaOutOfPlaneBend covalent force; also calls methods in the
# AmoebaHarmonicAngleGenerator to generate the AmoebaHarmonicAngleForce and
# AmoebaHarmonicInPlaneAngleForce
#=============================================================================================

class AmoebaOutOfPlaneBendGenerator:

    #=============================================================================================

    """An AmoebaOutOfPlaneBendGenerator constructs a AmoebaOutOfPlaneBendForce."""
    
    #=============================================================================================

    def __init__(self, forceField, type, cubic, quartic, pentic, sextic):

Peter Eastman's avatar
Peter Eastman committed
1656
1657
1658
1659
1660
1661
        self.forceField = forceField
        self.type = type
        self.cubic = cubic
        self.quartic = quartic
        self.pentic = pentic
        self.sextic = sextic
1662

Peter Eastman's avatar
Peter Eastman committed
1663
1664
1665
1666
        self.types1 = []
        self.types2 = []
        self.types3 = []
        self.types4 = []
1667

Peter Eastman's avatar
Peter Eastman committed
1668
1669
        self.ks = []
        self.hasBeenCalled = 0
1670
1671
1672
1673
1674
1675
1676
1677
1678
1679
1680
1681
1682
1683
1684
1685
1686
1687
1688
1689
1690
1691
1692
1693
1694
1695
1696
1697
1698
1699
1700
1701
1702
1703

    #=============================================================================================
    # Local version of findAtomTypes needed since class indices are 0 (i.e., not recognized)
    # for types3 and 4
    #=============================================================================================
    
    def findAtomTypes(self, forceField, node, num):
        """Parse the attributes on an XML tag to find the set of atom types for each atom it involves."""
        types = []
        attrib = node.attrib
        for i in range(num):
            if num == 1:
                suffix = ''
            else:
                suffix = str(i+1)
            classAttrib = 'class'+suffix
            if classAttrib in attrib:
                if attrib[classAttrib] in forceField._atomClasses:
                    types.append(forceField._atomClasses[attrib[classAttrib]])
                else:
                    types.append(set())
        return types

    #=============================================================================================

    @staticmethod
    def parseElement(element, forceField):

        #  <AmoebaOutOfPlaneBendForce type="ALLINGER" opbend-cubic="-0.014" opbend-quartic="5.6e-05" opbend-pentic="-7e-07" opbend-sextic="2.2e-08">
        #   <Angle class1="2" class2="1" class3="0" class4="0" k="0.0531474541591"/>
        #   <Angle class1="3" class2="1" class3="0" class4="0" k="0.0898536095496"/>
         
        # get global scalar parameters

Peter Eastman's avatar
Peter Eastman committed
1704
        generator = AmoebaOutOfPlaneBendGenerator(forceField, element.attrib['type'],
1705
1706
1707
                                                   float(element.attrib['opbend-cubic']),
                                                   float(element.attrib['opbend-quartic']),
                                                   float(element.attrib['opbend-pentic']),
Peter Eastman's avatar
Peter Eastman committed
1708
                                                   float(element.attrib['opbend-sextic']))
1709
1710
1711
1712

        forceField._forces.append(generator)

        for angle in element.findall('Angle'):
Peter Eastman's avatar
Peter Eastman committed
1713
            types = generator.findAtomTypes(forceField, angle, 4)
1714
1715
            if types is not None:

Peter Eastman's avatar
Peter Eastman committed
1716
1717
1718
1719
                generator.types1.append(types[0])
                generator.types2.append(types[1])
                generator.types3.append(types[2])
                generator.types4.append(types[3])
1720

Peter Eastman's avatar
Peter Eastman committed
1721
                generator.ks.append(float(angle.attrib['k']))
1722
1723
1724
1725
1726
1727

            else:
                outputString = "AmoebaOutOfPlaneBendGenerator: error getting types: %s %s %s %s." % (
                                    angle.attrib['class1'],
                                    angle.attrib['class2'],
                                    angle.attrib['class3'],
Peter Eastman's avatar
Peter Eastman committed
1728
1729
                                    angle.attrib['class4'])
                raise ValueError(outputString) 
1730
1731
1732
1733
1734
1735
1736
1737
    
    #=============================================================================================
    # Get middle atom in a angle
    # return index of middle atom or -1 if no middle is found
    # This method appears not to be needed since the angle[1] entry appears to always
    # be the middle atom. However, was unsure if this is guaranteed
    #=============================================================================================
    
Peter Eastman's avatar
Peter Eastman committed
1738
    def getMiddleAtom(self, angle, data):
1739
1740
1741

        # find atom shared by both bonds making up the angle

Peter Eastman's avatar
Peter Eastman committed
1742
        middleAtom = -1
1743
        for atomIndex in angle: 
Peter Eastman's avatar
Peter Eastman committed
1744
            isMiddle = 0
1745
1746
1747
            for bond in data.atomBonds[atomIndex]:
                atom1 = data.bonds[bond].atom1
                atom2 = data.bonds[bond].atom2
Peter Eastman's avatar
Peter Eastman committed
1748
                if (atom1 != atomIndex):
1749
1750
1751
                    partner = atom1
                else:
                    partner = atom2
Peter Eastman's avatar
Peter Eastman committed
1752
                if (partner == angle[0] or partner == angle[1] or partner == angle[2]): 
1753
1754
                    isMiddle += 1

Peter Eastman's avatar
Peter Eastman committed
1755
            if (isMiddle == 2):
1756
1757
1758
1759
1760
                return atomIndex
        return -1

    #=============================================================================================

Peter Eastman's avatar
Peter Eastman committed
1761
    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
1762

Peter Eastman's avatar
Peter Eastman committed
1763
1764
        verbose = 0
        if (self.hasBeenCalled):
1765
1766
1767
1768
1769
1770
1771
1772
1773
1774
1775
1776
1777
1778
1779
             return
        self.hasBeenCalled = 1

        # get force

        existing = [sys.getForce(i) for i in range(sys.getNumForces())]
        existing = [f for f in existing if type(f) == mm.AmoebaOutOfPlaneBendForce]
        if len(existing) == 0:
            force = mm.AmoebaOutOfPlaneBendForce()
            sys.addForce(force)
        else:
            force = existing[0]

        # set scalars

Peter Eastman's avatar
Peter Eastman committed
1780
1781
1782
1783
        force.setAmoebaGlobalOutOfPlaneBendCubic(  self.cubic)
        force.setAmoebaGlobalOutOfPlaneBendQuartic(self.quartic)
        force.setAmoebaGlobalOutOfPlaneBendPentic( self.pentic)
        force.setAmoebaGlobalOutOfPlaneBendSextic( self.sextic)
1784

Peter Eastman's avatar
Peter Eastman committed
1785
1786
        count = 0
        opBondCount = 0
1787
1788
1789
1790

        # this hash is used to insure the out-of-plane-bend bonds
        # are only added once

Peter Eastman's avatar
Peter Eastman committed
1791
        skipAtoms = dict()
1792
1793
1794
1795

        # these lists are used in the partitioning of the angles into
        # angle and inPlane angles

Peter Eastman's avatar
Peter Eastman committed
1796
1797
        inPlaneAngles = []
        nonInPlaneAngles = []
1798
        nonInPlaneAnglesConstrained = []
Peter Eastman's avatar
Peter Eastman committed
1799
        idealAngles = []*len(data.angles)
1800
1801
1802

        for (angle, isConstrained) in zip(data.angles, data.isAngleConstrained):

Peter Eastman's avatar
Peter Eastman committed
1803
1804
1805
            middleAtom = self.getMiddleAtom(angle, data)
            if (middleAtom > -1):
                middleType = data.atomType[data.atoms[middleAtom]]
1806
1807
                middleCovalency = len(data.atomBonds[middleAtom])
            else:
Peter Eastman's avatar
Peter Eastman committed
1808
                middleType = -1
1809
1810
1811
1812
1813
1814
1815
1816
                middleCovalency = -1

            # if middle atom has covalency of 3 and 
            # the types of the middle atom and the partner atom (atom bonded to
            # middle atom, but not in angle) match types2 and types2, then
            # three out-of-plane bend angles are generated. Three in-plane angle 
            # are also generated. If the conditions are not satisfied the angle is generic angle (not a in-plane angle)

Peter Eastman's avatar
Peter Eastman committed
1817
            if (middleAtom > -1 and middleCovalency == 3 and middleAtom not in skipAtoms):
1818

Peter Eastman's avatar
Peter Eastman committed
1819
1820
1821
1822
                partners = []
                partnerSet = set()
                partnerTypes = []
                partnerK = []
1823
1824
1825
1826

                for bond in data.atomBonds[middleAtom]:
                    atom1 = data.bonds[bond].atom1
                    atom2 = data.bonds[bond].atom2
Peter Eastman's avatar
Peter Eastman committed
1827
                    if (atom1 != middleAtom):
1828
1829
1830
1831
1832
1833
1834
1835
                        partner = atom1
                    else:
                        partner = atom2

                    partnerType = data.atomType[data.atoms[partner]]
                    for i in range(len(self.types1)):
                        types1 = self.types1[i]
                        types2 = self.types2[i]
Peter Eastman's avatar
Peter Eastman committed
1836
1837
1838
1839
1840
                        if (middleType in types2 and partnerType in types1):
                            partners.append(partner)
                            partnerSet.add(partner)
                            partnerTypes.append(partnerType)
                            partnerK.append(self.ks[i])
1841
             
Peter Eastman's avatar
Peter Eastman committed
1842
                if (len(partners) == 3):
1843
1844

                    opBondCount += 3
Peter Eastman's avatar
Peter Eastman committed
1845
                    if (verbose):
1846
1847
                        print "%5d Opbend: %d type=%s cov=%d " % (opBondCount, middleAtom, middleType, middleCovalency)

Peter Eastman's avatar
Peter Eastman committed
1848
1849
1850
                    force.addOutOfPlaneBend(partners[0], middleAtom, partners[1], partners[2], partnerK[2])
                    force.addOutOfPlaneBend(partners[0], middleAtom, partners[2], partners[1], partnerK[1])
                    force.addOutOfPlaneBend(partners[1], middleAtom, partners[2], partners[0], partnerK[0])
1851
1852

                    skipAtoms[middleAtom] = set()
Peter Eastman's avatar
Peter Eastman committed
1853
1854
1855
1856
1857
1858
1859
1860
1861
                    skipAtoms[middleAtom].add(partners[0])
                    skipAtoms[middleAtom].add(partners[1])
                    skipAtoms[middleAtom].add(partners[2])

                    angleDict = {}
                    angleList = []
                    angleList.append(angle[0])
                    angleList.append(angle[1])
                    angleList.append(angle[2])
1862
1863
                    angleDict['angle'] = angleList

Peter Eastman's avatar
Peter Eastman committed
1864
                    angleDict['isConstrained'] = 0
1865

Peter Eastman's avatar
Peter Eastman committed
1866
1867
1868
1869
                    angleSet = set()
                    angleSet.add(angle[0])
                    angleSet.add(angle[1])
                    angleSet.add(angle[2])
1870
1871

                    for atomIndex in partnerSet:
Peter Eastman's avatar
Peter Eastman committed
1872
1873
                        if (atomIndex not in angleSet):
                            angleList.append(atomIndex)
1874

Peter Eastman's avatar
Peter Eastman committed
1875
1876
                    if (verbose):
                        print "%5d Opbend: middle=%d inPlane1 %d %s" % (opBondCount, middleAtom, len(angleList), str(angleList))
1877

Peter Eastman's avatar
Peter Eastman committed
1878
                    inPlaneAngles.append(angleDict)
1879
1880

                else:
Peter Eastman's avatar
Peter Eastman committed
1881
1882
1883
1884
1885
1886
1887
                    if (verbose):
                        print "%5d Opbend: %d type=%s cov=%d  xxx" % (opBondCount, middleAtom, middleType, middleCovalency)
                        print "%5d Opbend: middle=%d inPlane1 %d %s" % (opBondCount, middleAtom, len(angleList), str(angleList))
                    angleDict = {}
                    angleDict['angle'] = angle
                    angleDict['isConstrained'] = isConstrained
                    nonInPlaneAngles.append(angleDict)
1888
            else:
Peter Eastman's avatar
Peter Eastman committed
1889
                if (middleAtom > -1 and middleCovalency == 3 and middleAtom in skipAtoms):
1890

Peter Eastman's avatar
Peter Eastman committed
1891
                    partnerSet = skipAtoms[middleAtom]
1892
                  
Peter Eastman's avatar
Peter Eastman committed
1893
                    angleDict = {}
1894

Peter Eastman's avatar
Peter Eastman committed
1895
1896
1897
1898
1899
                    angleList = []
                    angleList.append(angle[0])
                    angleList.append(angle[1])
                    angleList.append(angle[2])
                    angleDict['angle'] = angleList
1900

Peter Eastman's avatar
Peter Eastman committed
1901
                    angleDict['isConstrained'] = isConstrained
1902

Peter Eastman's avatar
Peter Eastman committed
1903
1904
1905
1906
                    angleSet = set()
                    angleSet.add(angle[0])
                    angleSet.add(angle[1])
                    angleSet.add(angle[2])
1907
1908

                    for atomIndex in partnerSet:
Peter Eastman's avatar
Peter Eastman committed
1909
1910
                        if (atomIndex not in angleSet):
                            angleList.append(atomIndex)
1911

Peter Eastman's avatar
Peter Eastman committed
1912
1913
                    if (verbose):
                        print "%5d Opbend: middle=%d inPlane2 %d angleList=%s partnerSet=%s" % (opBondCount, middleAtom, len(angleList), str(angleList), str(partnerSet))
1914

Peter Eastman's avatar
Peter Eastman committed
1915
                    inPlaneAngles.append(angleDict)
1916
1917

                else:
Peter Eastman's avatar
Peter Eastman committed
1918
1919
                    angleDict = {}
                    angleDict['angle'] = angle
1920
                    angleDict['isConstrained'] = isConstrained
Peter Eastman's avatar
Peter Eastman committed
1921
                    nonInPlaneAngles.append(angleDict)
1922
1923
1924
1925
1926
1927

            count += 1

        # get AmoebaHarmonicAngleGenerator and add AmoebaHarmonicAngle and AmoebaHarmonicInPlaneAngle forces

        for force in self.forceField._forces:
Peter Eastman's avatar
Peter Eastman committed
1928
1929
1930
1931
1932
            if (force.__class__.__name__ == 'AmoebaHarmonicAngleGenerator'): 
                force.createForcePostOpBendAngle(sys, data, nonbondedMethod, nonbondedCutoff, nonInPlaneAngles, args)
                force.createForcePostOpBendInPlaneAngle(sys, data, nonbondedMethod, nonbondedCutoff, inPlaneAngles, args)
            if (force.__class__.__name__ == 'AmoebaHarmonicBondGenerator'): 
                force.createForce(sys, data, nonbondedMethod, nonbondedCutoff, args)
1933
1934

        for force in self.forceField._forces:
Peter Eastman's avatar
Peter Eastman committed
1935
            if (force.__class__.__name__ == 'AmoebaStretchBendGenerator'): 
1936
                for angleDict in inPlaneAngles:
Peter Eastman's avatar
Peter Eastman committed
1937
1938
                    nonInPlaneAngles.append(angleDict)
                force.createForcePostAmoebaHarmonicBondForce(sys, data, nonbondedMethod, nonbondedCutoff, nonInPlaneAngles, args)
1939
1940
1941
1942
1943
1944
1945
1946
1947
1948
1949
1950
1951
1952
1953

parsers["AmoebaOutOfPlaneBendForce"] = AmoebaOutOfPlaneBendGenerator.parseElement

#=============================================================================================

class AmoebaTorsionGenerator:

    #=============================================================================================

    """An AmoebaTorsionGenerator constructs a AmoebaTorsionForce."""

    #=============================================================================================

    def __init__(self, torsionUnit):

Peter Eastman's avatar
Peter Eastman committed
1954
        self.torsionUnit = torsionUnit
1955

Peter Eastman's avatar
Peter Eastman committed
1956
1957
1958
1959
        self.types1 = []
        self.types2 = []
        self.types3 = []
        self.types4 = []
1960

Peter Eastman's avatar
Peter Eastman committed
1961
1962
1963
        self.t1 = []
        self.t2 = []
        self.t3 = []
1964
1965
1966
1967
1968
1969
1970
1971
1972
1973
1974
        self.hasBeenCalled = 0
    
    #=============================================================================================

    @staticmethod
    def parseElement(element, forceField):

        #  <AmoebaTorsionForce torsionUnit="0.5">
        #   <Torsion class1="3" class2="1" class3="2" class4="3"   amp1="0.0" angle1="0.0"   amp2="0.0" angle2="3.14159265359"   amp3="0.0" angle3="0.0" />
        #   <Torsion class1="3" class2="1" class3="2" class4="6"   amp1="0.0" angle1="0.0"   amp2="0.0" angle2="3.14159265359"   amp3="-0.263592" angle3="0.0" />
         
Peter Eastman's avatar
Peter Eastman committed
1975
        generator = AmoebaTorsionGenerator(float(element.attrib['torsionUnit']))
1976
1977
1978
1979
1980
1981
1982
1983
1984
1985
1986
1987
1988
1989
1990
        forceField._forces.append(generator)

        # collect particle classes and t1,t2,t3,
        # where ti=[amplitude_i,angle_i]

        for torsion in element.findall('Torsion'):
            types = forceField._findAtomTypes(torsion, 4)
            if types is not None:

                generator.types1.append(types[0])
                generator.types2.append(types[1])
                generator.types3.append(types[2])
                generator.types4.append(types[3])

                for ii in range(1,4):
Peter Eastman's avatar
Peter Eastman committed
1991
1992
                    tInfo = []
                    suffix = str(ii)
1993
1994
1995
1996
1997
1998
                    ampName = 'amp' + suffix
                    tInfo.append(float(torsion.attrib[ampName]))

                    angName = 'angle' + suffix
                    tInfo.append(float(torsion.attrib[angName]))

Peter Eastman's avatar
Peter Eastman committed
1999
2000
2001
2002
2003
2004
                    if (ii == 1):
                        generator.t1.append(tInfo)
                    elif (ii == 2):
                        generator.t2.append(tInfo)
                    elif (ii == 3):
                        generator.t3.append(tInfo)
2005
2006
2007
2008
2009
2010

            else:
                outputString = "AmoebaTorsionGenerator: error getting types: %s %s %s %s" % (
                                    stretchBend.attrib['class1'],
                                    stretchBend.attrib['class2'],
                                    stretchBend.attrib['class3'],
Peter Eastman's avatar
Peter Eastman committed
2011
2012
                                    stretchBend.attrib['class4'])
                raise ValueError(outputString) 
2013
2014
2015
2016
2017
    
    #=============================================================================================

    def createForce(self, sys, data, nontorsionedMethod, nontorsionedCutoff, args):

Peter Eastman's avatar
Peter Eastman committed
2018
2019
        verbose = 0
        if (self.hasBeenCalled):
2020
2021
2022
2023
2024
2025
2026
2027
2028
2029
             return
        self.hasBeenCalled = 1

        existing = [sys.getForce(i) for i in range(sys.getNumForces())]
        existing = [f for f in existing if type(f) == mm.AmoebaTorsionForce]
        if len(existing) == 0:
            force = mm.AmoebaTorsionForce()
            sys.addForce(force)
        else:
            force = existing[0]
Peter Eastman's avatar
Peter Eastman committed
2030
        if (verbose):
2031
2032
2033
2034
2035
2036
2037
2038
2039
            print "In AmoebaTorsionGenerator torsions=%d " % (len(data.propers))
        count = 0
        for torsion in data.propers:

            type1 = data.atomType[data.atoms[torsion[0]]]
            type2 = data.atomType[data.atoms[torsion[1]]]
            type3 = data.atomType[data.atoms[torsion[2]]]
            type4 = data.atomType[data.atoms[torsion[3]]]

Peter Eastman's avatar
Peter Eastman committed
2040
            hit = 0
2041
2042
2043
2044
2045
2046
2047
2048
2049
            for i in range(len(self.types1)):

                types1 = self.types1[i]
                types2 = self.types2[i]
                types3 = self.types3[i]
                types4 = self.types4[i]

                # match types in forward or reverse direction

Peter Eastman's avatar
Peter Eastman committed
2050
2051
                if (type1 in types1 and type2 in types2 and type3 in types3 and type4 in types4) or (type4 in types1 and type3 in types2 and type2 in types3 and type1 in types4):
                    if (verbose):
2052
2053
2054
                        print "AmoebaTorsionGenerator %5d [%5d %5d %5d %5d] [%5s %5s %5s %5s]" % (
                              count, torsion[0], torsion[1], torsion[2], torsion[3],
                              type1, type2, type3, type4)
Peter Eastman's avatar
Peter Eastman committed
2055
2056
                    hit = 1
                    force.addTorsion(torsion[0], torsion[1], torsion[2], torsion[3], self.t1[i],  self.t2[i], self.t3[i])
2057
                    break
Peter Eastman's avatar
Peter Eastman committed
2058
            if (hit == 0): 
2059
2060
2061
2062
                print "AmoebaTorsionGenerator missing: %5d %s %s " % (count, type1, type2, type3, type4)

            count += 1

Peter Eastman's avatar
Peter Eastman committed
2063
        if (verbose):
2064
2065
2066
2067
2068
2069
2070
2071
2072
2073
2074
2075
2076
2077
2078
2079
            print "AmoebaTorsionGenerator number of torsions added=%d " % (force.getNumTorsions())

parsers["AmoebaTorsionForce"] = AmoebaTorsionGenerator.parseElement

#=============================================================================================

class AmoebaPiTorsionGenerator:

    #=============================================================================================

    """An AmoebaPiTorsionGenerator constructs a AmoebaPiTorsionForce."""

    #=============================================================================================
    
    def __init__(self, piTorsionUnit):
        self.piTorsionUnit = piTorsionUnit 
Peter Eastman's avatar
Peter Eastman committed
2080
2081
2082
        self.types1 = []
        self.types2 = []
        self.k = []
2083
2084
2085
2086
2087
2088
2089
2090
2091
2092
        self.hasBeenCalled = 0
    
    #=============================================================================================

    @staticmethod
    def parseElement(element, forceField):

        #  <AmoebaPiTorsionForce piTorsionUnit="1.0">
        #   <PiTorsion class1="1" class2="3" k="28.6604" />

Peter Eastman's avatar
Peter Eastman committed
2093
        generator = AmoebaPiTorsionGenerator(float(element.attrib['piTorsionUnit']))
2094
2095
2096
2097
2098
2099
2100
2101
2102
2103
2104
        forceField._forces.append(generator)

        for piTorsion in element.findall('PiTorsion'):
            types = forceField._findAtomTypes(piTorsion, 2)
            if types is not None:
                generator.types1.append(types[0])
                generator.types2.append(types[1])
                generator.k.append(float(piTorsion.attrib['k']))
            else:
                outputString = "AmoebaPiTorsionGenerator: error getting types: %s %s " % (
                                    piTorsion.attrib['class1'],
Peter Eastman's avatar
Peter Eastman committed
2105
2106
                                    piTorsion.attrib['class2'])
                raise ValueError(outputString) 
2107
2108
2109
2110
2111
    
    #=============================================================================================

    def createForce(self, sys, data, nonpiTorsionedMethod, nonpiTorsionedCutoff, args):

Peter Eastman's avatar
Peter Eastman committed
2112
2113
        verbose = 0
        if (self.hasBeenCalled):
2114
2115
2116
2117
2118
2119
2120
2121
2122
2123
2124
2125
2126
2127
2128
2129
2130
2131
2132
2133
             return
        self.hasBeenCalled = 1

        existing = [sys.getForce(i) for i in range(sys.getNumForces())]
        existing = [f for f in existing if type(f) == mm.AmoebaPiTorsionForce]

        if len(existing) == 0:
            force = mm.AmoebaPiTorsionForce()
            sys.addForce(force)
        else:
            force = existing[0]

        count = 0
        for bond in data.bonds:

            # search for bonds with both atoms in bond having covalency == 3
 
            atom1 = bond.atom1
            atom2 = bond.atom2
 
Peter Eastman's avatar
Peter Eastman committed
2134
            if (len(data.atomBonds[atom1]) == 3 and len(data.atomBonds[atom1]) == 3):
2135
2136
2137
2138
2139
2140
2141
2142
2143
2144

                type1 = data.atomType[data.atoms[atom1]]
                type2 = data.atomType[data.atoms[atom2]]

                for i in range(len(self.types1)):

                   types1 = self.types1[i]
                   types2 = self.types2[i]

                   if (type1 in types1 and type2 in types2) or (type1 in types2 and type2 in types1):
Peter Eastman's avatar
Peter Eastman committed
2145
                       if (verbose):
2146
2147
2148
2149
2150
2151
2152
2153
2154
2155
2156
2157
2158
2159
2160
2161
                            print "AmoebaPiTorsionGenerator %5d %5d %5d [%5s %5s] %15.6f" % (count, atom1, atom2, type1, type2, self.k[i])

                       # piTorsionAtom1, piTorsionAtom2 are the atoms bonded to atom1, excluding atom2 
                       # piTorsionAtom5, piTorsionAtom6 are the atoms bonded to atom2, excluding atom1 

                       piTorsionAtom1 = -1
                       piTorsionAtom2 = -1
                       piTorsionAtom3 = atom1

                       piTorsionAtom4 = atom2
                       piTorsionAtom5 = -1
                       piTorsionAtom6 = -1

                       for bond in data.atomBonds[atom1]:
                           bondedAtom1 = data.bonds[bond].atom1
                           bondedAtom2 = data.bonds[bond].atom2
Peter Eastman's avatar
Peter Eastman committed
2162
                           if (bondedAtom1 != atom1):
2163
2164
2165
                               b1 = bondedAtom1
                           else:
                               b1 = bondedAtom2
Peter Eastman's avatar
Peter Eastman committed
2166
2167
                           if (b1 != atom2):
                               if (piTorsionAtom1 == -1):
2168
2169
2170
2171
2172
2173
2174
                                   piTorsionAtom1 = b1 
                               else:
                                   piTorsionAtom2 = b1

                       for bond in data.atomBonds[atom2]:
                           bondedAtom1 = data.bonds[bond].atom1
                           bondedAtom2 = data.bonds[bond].atom2
Peter Eastman's avatar
Peter Eastman committed
2175
                           if (bondedAtom1 != atom2):
2176
2177
2178
2179
                               b1 = bondedAtom1
                           else:
                               b1 = bondedAtom2

Peter Eastman's avatar
Peter Eastman committed
2180
2181
                           if (b1 != atom1):
                               if (piTorsionAtom5 == -1):
2182
2183
2184
2185
                                   piTorsionAtom5 = b1 
                               else:
                                   piTorsionAtom6 = b1
    
Peter Eastman's avatar
Peter Eastman committed
2186
                       force.addPiTorsion(piTorsionAtom1, piTorsionAtom2, piTorsionAtom3, piTorsionAtom4, piTorsionAtom5, piTorsionAtom6, self.k[i])
2187
2188

            count += 1
Peter Eastman's avatar
Peter Eastman committed
2189
        if (verbose):
2190
2191
2192
2193
2194
2195
2196
2197
2198
2199
2200
2201
2202
2203
            print "AmoebaPiTorsionGenerator number of pi-torsions added=%d " % (force.getNumPiTorsions())

parsers["AmoebaPiTorsionForce"] = AmoebaPiTorsionGenerator.parseElement

#=============================================================================================

class AmoebaTorsionTorsionGenerator:

    #=============================================================================================

    """An AmoebaTorsionTorsionGenerator constructs a AmoebaTorsionTorsionForce."""
    
    #=============================================================================================

Peter Eastman's avatar
Peter Eastman committed
2204
    def __init__(self):
2205

Peter Eastman's avatar
Peter Eastman committed
2206
2207
2208
2209
2210
        self.types1 = []
        self.types2 = []
        self.types3 = []
        self.types4 = []
        self.types5 = []
2211

Peter Eastman's avatar
Peter Eastman committed
2212
        self.gridIndex = []
2213

Peter Eastman's avatar
Peter Eastman committed
2214
        self.grids = []
2215
2216
2217
2218
2219
2220
2221
2222

        self.hasBeenCalled = 0
    
    #=============================================================================================

    @staticmethod
    def parseElement(element, forceField):

Peter Eastman's avatar
Peter Eastman committed
2223
        generator = AmoebaTorsionTorsionGenerator()
2224
2225
2226
2227
2228
2229
2230
        forceField._forces.append(generator)
        maxGridIndex = -1

        # <AmoebaTorsionTorsionForce >
        # <TorsionTorsion class1="3" class2="1" class3="2" class4="3" class5="1" grid="0" nx="25" ny="25" />

        for torsionTorsion in element.findall('TorsionTorsion'):
Peter Eastman's avatar
Peter Eastman committed
2231
            types = forceField._findAtomTypes(torsionTorsion, 5)
2232
2233
2234
2235
2236
2237
2238
2239
            if types is not None:

                generator.types1.append(types[0])
                generator.types2.append(types[1])
                generator.types3.append(types[2])
                generator.types4.append(types[3])
                generator.types5.append(types[4])

Peter Eastman's avatar
Peter Eastman committed
2240
2241
                gridIndex = int(torsionTorsion.attrib['grid'])
                if (gridIndex > maxGridIndex):
2242
2243
2244
2245
2246
2247
2248
2249
2250
                    maxGridIndex = gridIndex

                generator.gridIndex.append(gridIndex)
            else:
                outputString = "AmoebaTorsionTorsionGenerator: error getting types: %s %s %s %s %s" % (
                                    torsionTorsion.attrib['class1'],
                                    torsionTorsion.attrib['class2'],
                                    torsionTorsion.attrib['class3'],
                                    torsionTorsion.attrib['class4'],
Peter Eastman's avatar
Peter Eastman committed
2251
2252
                                    torsionTorsion.attrib['class5'] )
                raise ValueError(outputString) 
2253
2254
2255
2256
2257
2258
2259
2260
2261
2262
2263
2264
2265
2266
2267
2268
2269
2270
2271
    
        # load grid

        # xml source

        # <TorsionTorsionGrid grid="0" nx="25" ny="25" >
        # <Grid angle1="-180.00" angle2="-180.00" f="0.0" fx="2.31064374824e-05" fy="0.0" fxy="-0.0052801799672" />
        # <Grid angle1="-165.00" angle2="-180.00" f="-0.66600912" fx="-0.06983370052" fy="-0.075058725744" fxy="-0.0044462732032" />

        # output grid:

        #     grid[x][y][0] = x value 
        #     grid[x][y][1] = y value 
        #     grid[x][y][2] = function value 
        #     grid[x][y][3] = dfdx value 
        #     grid[x][y][4] = dfdy value 
        #     grid[x][y][5] = dfd(xy) value 

        maxGridIndex    += 1
Peter Eastman's avatar
Peter Eastman committed
2272
        generator.grids = maxGridIndex*[]
2273
2274
        for torsionTorsionGrid in element.findall('TorsionTorsionGrid'):

Peter Eastman's avatar
Peter Eastman committed
2275
2276
2277
            gridIndex = int(torsionTorsionGrid.attrib[ "grid"])
            nx = int(torsionTorsionGrid.attrib[ "nx"])
            ny = int(torsionTorsionGrid.attrib[ "ny"])
2278

Peter Eastman's avatar
Peter Eastman committed
2279
2280
            grid = []
            gridCol = []
2281
2282
2283
2284
2285

            gridColIndex = 0

            for gridEntry in torsionTorsionGrid.findall('Grid'):

Peter Eastman's avatar
Peter Eastman committed
2286
2287
2288
2289
2290
2291
2292
2293
                gridRow = []
                gridRow.append(float(gridEntry.attrib['angle1']))
                gridRow.append(float(gridEntry.attrib['angle2']))
                gridRow.append(float(gridEntry.attrib['f']))
                gridRow.append(float(gridEntry.attrib['fx']))
                gridRow.append(float(gridEntry.attrib['fy']))
                gridRow.append(float(gridEntry.attrib['fxy']))
                gridCol.append(gridRow)
2294
2295

                gridColIndex  += 1
Peter Eastman's avatar
Peter Eastman committed
2296
2297
2298
                if (gridColIndex == nx):
                    grid.append(gridCol)
                    gridCol = []
2299
2300
2301
                    gridColIndex = 0

            
Peter Eastman's avatar
Peter Eastman committed
2302
2303
            if (gridIndex == len(generator.grids)):
                generator.grids.append(grid)
2304
            else:
Peter Eastman's avatar
Peter Eastman committed
2305
2306
                while(len(generator.grids) < gridIndex):
                    generator.grids.append([])
2307
2308
2309
2310
                generator.grids[gridIndex] = grid

    #=============================================================================================

Peter Eastman's avatar
Peter Eastman committed
2311
    def getChiralAtomIndex(self, data, sys, atomB, atomC, atomD):
2312
2313
2314
2315
2316
2317
2318
2319

        chiralAtomIndex = -1

        # if atomC has four bonds, find the
        # two bonds that do not include atomB and atomD
        # set chiralAtomIndex to one of these, if they are
        # not the same atom(type/mass)

Peter Eastman's avatar
Peter Eastman committed
2320
        if (len(data.atomBonds[atomC]) == 4):
2321
2322
2323
2324
2325
            atomE = -1
            atomF = -1
            for bond in data.atomBonds[atomC]:
                bondedAtom1 = data.bonds[bond].atom1
                bondedAtom2 = data.bonds[bond].atom2
Peter Eastman's avatar
Peter Eastman committed
2326
2327
                hit = -1
                if (  bondedAtom1 == atomC and bondedAtom2 != atomB and bondedAtom2 != atomD):
2328
                    hit = bondedAtom2
Peter Eastman's avatar
Peter Eastman committed
2329
                elif (bondedAtom2 == atomC and bondedAtom1 != atomB and bondedAtom1 != atomD):
2330
2331
                    hit = bondedAtom1

Peter Eastman's avatar
Peter Eastman committed
2332
2333
                if (hit > -1):
                    if (atomE == -1):
2334
2335
2336
2337
2338
2339
                        atomE = hit
                    else:
                        atomF = hit
       
            # raise error if atoms E or F not found

Peter Eastman's avatar
Peter Eastman committed
2340
2341
2342
            if (atomE == -1 or atomF == -1):
                outputString = "getChiralAtomIndex: error getting bonded partners of atomC=%s %d %s" % (atomC.name, atomC.resiude.index, atomC.resiude.name,)
                raise ValueError(outputString) 
2343
2344
2345
2346
2347

            # check for different type/mass between atoms E & F

            typeE = int(data.atomType[data.atoms[atomE]])
            typeF = int(data.atomType[data.atoms[atomF]])
Peter Eastman's avatar
Peter Eastman committed
2348
            if (typeE > typeF):
2349
                chiralAtomIndex = atomE 
Peter Eastman's avatar
Peter Eastman committed
2350
            if (typeF > typeE):
2351
2352
                chiralAtomIndex = atomF 

Peter Eastman's avatar
Peter Eastman committed
2353
2354
2355
            massE = sys.getParticleMass(atomE)/unit.dalton
            massF = sys.getParticleMass(atomE)/unit.dalton
            if (massE > massF):
2356
                chiralAtomIndex = massE 
Peter Eastman's avatar
Peter Eastman committed
2357
            if (massF > massE):
2358
2359
2360
2361
2362
2363
2364
2365
                chiralAtomIndex = massF 

        return chiralAtomIndex

    #=============================================================================================
 
    def createForce(self, sys, data, nonpiTorsionedMethod, nonpiTorsionedCutoff, args):

Peter Eastman's avatar
Peter Eastman committed
2366
2367
        verbose = 0
        if (self.hasBeenCalled):
2368
2369
2370
2371
2372
2373
2374
2375
2376
2377
2378
2379
2380
2381
2382
2383
2384
2385
2386
2387
2388
2389
2390
2391
             return
        self.hasBeenCalled = 1

        existing = [sys.getForce(i) for i in range(sys.getNumForces())]
        existing = [f for f in existing if type(f) == mm.AmoebaTorsionTorsionForce]

        if len(existing) == 0:
            force = mm.AmoebaTorsionTorsionForce()
            sys.addForce(force)
        else:
            force = existing[0]

        count = 0
        for angle in data.angles:

            # search for bitorsions; based on TINKER subroutine bitors() 
 
            ib = angle[0]
            ic = angle[1]
            id = angle[2]

            for bondIndex in data.atomBonds[ib]:
                bondedAtom1 = data.bonds[bondIndex].atom1
                bondedAtom2 = data.bonds[bondIndex].atom2
Peter Eastman's avatar
Peter Eastman committed
2392
                if (bondedAtom1 != ib):
2393
2394
2395
2396
                    ia = bondedAtom1
                else:
                    ia = bondedAtom2

Peter Eastman's avatar
Peter Eastman committed
2397
                if (ia != ic and ia != id):
2398
2399
2400
                    for bondIndex in data.atomBonds[id]:
                        bondedAtom1 = data.bonds[bondIndex].atom1
                        bondedAtom2 = data.bonds[bondIndex].atom2
Peter Eastman's avatar
Peter Eastman committed
2401
                        if (bondedAtom1 != id):
2402
2403
2404
2405
                            ie = bondedAtom1
                        else:
                            ie = bondedAtom2

Peter Eastman's avatar
Peter Eastman committed
2406
                        if (ie != ic and ie != ib and ie != ia):
2407
2408
2409
2410
2411
2412
2413
2414
2415
2416
2417
2418
2419
2420
2421
2422
2423
2424
2425
2426

                            # found candidate set of atoms
                            # check if types match in order or reverse order

                            type1 = data.atomType[data.atoms[ia]]
                            type2 = data.atomType[data.atoms[ib]]
                            type3 = data.atomType[data.atoms[ic]]
                            type4 = data.atomType[data.atoms[id]]
                            type5 = data.atomType[data.atoms[ie]]

                            for i in range(len(self.types1)):

                                types1 = self.types1[i]
                                types2 = self.types2[i]
                                types3 = self.types3[i]
                                types4 = self.types4[i]
                                types5 = self.types5[i]

                                # match in order

Peter Eastman's avatar
Peter Eastman committed
2427
2428
2429
                                if (type1 in types1 and type2 in types2 and type3 in types3 and type4 in types4 and type5 in types5):
                                    chiralAtomIndex = self.getChiralAtomIndex(data, sys, ib, ic, id)
                                    force.addTorsionTorsion(ia, ib, ic, id, ie, chiralAtomIndex, self.gridIndex[i])
2430
2431
2432

                                # match in reverse order

Peter Eastman's avatar
Peter Eastman committed
2433
2434
2435
                                if (type5 in types1 and type4 in types2 and type3 in types3 and type2 in types4 and type1 in types5):
                                    chiralAtomIndex = self.getChiralAtomIndex(data, sys, ib, ic, id)
                                    force.addTorsionTorsion(ie, id, ic, ib, ia, chiralAtomIndex, self.gridIndex[i])
2436
2437
2438
2439

        # set grids

        for (index, grid) in enumerate(self.grids):
Peter Eastman's avatar
Peter Eastman committed
2440
            force.setTorsionTorsionGrid(index, grid)
2441
 
Peter Eastman's avatar
Peter Eastman committed
2442
        if (verbose):
2443
2444
2445
2446
2447
2448
2449
2450
2451
2452
2453
2454
2455
2456
            print "AmoebaTorsionTorsionGenerator number of bitorsions added=%d " % (force.getNumTorsionTorsions())
            print "AmoebaTorsionTorsionGenerator number of grids      added=%d " % (force.getNumTorsionTorsionGrids())

parsers["AmoebaTorsionTorsionForce"] = AmoebaTorsionTorsionGenerator.parseElement

#=============================================================================================

class AmoebaStretchBendGenerator:
    """An AmoebaStretchBendGenerator constructs a AmoebaStretchBendForce."""
    
    #=============================================================================================

    def __init__(self):

Peter Eastman's avatar
Peter Eastman committed
2457
2458
2459
        self.types1 = []
        self.types2 = []
        self.types3 = []
2460

Peter Eastman's avatar
Peter Eastman committed
2461
2462
        self.k1 = []
        self.k2 = []
2463
2464
2465
2466
2467
2468
        self.hasBeenCalled = 0
    
    #=============================================================================================

    @staticmethod
    def parseElement(element, forceField):
Peter Eastman's avatar
Peter Eastman committed
2469
        generator = AmoebaStretchBendGenerator()
2470
2471
2472
2473
2474
2475
2476
2477
2478
2479
2480
2481
2482
2483
2484
2485
2486
2487
2488
2489
2490
        forceField._forces.append(generator)

        # <AmoebaStretchBendForce stretchBendUnit="1.0">
        # <StretchBend class1="2" class2="1" class3="3" k1="5.25776946506" k2="5.25776946506" />
        # <StretchBend class1="2" class2="1" class3="4" k1="3.14005676385" k2="3.14005676385" />

        for stretchBend in element.findall('StretchBend'):
            types = forceField._findAtomTypes(stretchBend, 3)
            if types is not None:

                generator.types1.append(types[0])
                generator.types2.append(types[1])
                generator.types3.append(types[2])

                generator.k1.append(float(stretchBend.attrib['k1']))
                generator.k2.append(float(stretchBend.attrib['k2']))

            else:
                outputString = "AmoebaStretchBendGenerator : error getting types: %s %s %s" % (
                                    stretchBend.attrib['class1'],
                                    stretchBend.attrib['class2'],
Peter Eastman's avatar
Peter Eastman committed
2491
2492
                                    stretchBend.attrib['class3'])
                raise ValueError(outputString) 
2493
2494
2495
2496
2497
2498
2499
2500
2501
2502
2503
2504
    
    #=============================================================================================

    # The setup of this force is dependent on AmoebaHarmonicBondForce and AmoebaHarmonicAngleForce 
    # having been called since the ideal bond lengths and angle are needed here.
    # As a conseqeunce, createForce() is not implemented since it is not guaranteed that the generator for
    # AmoebaHarmonicBondForce and AmoebaHarmonicAngleForce have been called prior to AmoebaStretchBendGenerator(). 
    # Instead, createForcePostAmoebaHarmonicBondForce() is called 
    # after the generators for AmoebaHarmonicBondForce and AmoebaHarmonicAngleForce have been called

    #=============================================================================================

Peter Eastman's avatar
Peter Eastman committed
2505
    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
2506
2507
2508
2509
2510
2511
2512
2513
2514
2515
        pass

    #=============================================================================================

    # Note: request for constrained bonds is ignored.

    #=============================================================================================

    def createForcePostAmoebaHarmonicBondForce(self, sys, data, nonbondedMethod, nonbondedCutoff, angleList, args):

Peter Eastman's avatar
Peter Eastman committed
2516
2517
        verbose = 0
        if (self.hasBeenCalled):
2518
             return
Peter Eastman's avatar
Peter Eastman committed
2519
        self.hasBeenCalled = 1
2520
2521
2522
2523
2524
2525
2526
2527
2528

        existing = [sys.getForce(i) for i in range(sys.getNumForces())]
        existing = [f for f in existing if type(f) == mm.AmoebaStretchBendForce]
        if len(existing) == 0:
            force = mm.AmoebaStretchBendForce()
            sys.addForce(force)
        else:
            force = existing[0]

Peter Eastman's avatar
Peter Eastman committed
2529
        if (verbose):
2530
2531
            print "In AmoebaStretchBendGenerator bonds=%d " % (len(data.bonds))

Peter Eastman's avatar
Peter Eastman committed
2532
        count = 0
2533
2534
2535
        for angleDict in angleList:

            angle = angleDict['angle']
Peter Eastman's avatar
Peter Eastman committed
2536
            if ('isConstrained' in angleDict):
2537
2538
2539
2540
2541
2542
2543
2544
                isConstrained = angleDict['isConstrained']
            else:
                isConstrained = 0

            type1 = data.atomType[data.atoms[angle[0]]]
            type2 = data.atomType[data.atoms[angle[1]]]
            type3 = data.atomType[data.atoms[angle[2]]]

Peter Eastman's avatar
Peter Eastman committed
2545
2546
            hit = 0
            radian = 57.2957795130
2547
2548
2549
2550
2551
2552
2553
2554
2555
2556
            for i in range(len(self.types1)):

                types1 = self.types1[i]
                types2 = self.types2[i]
                types3 = self.types3[i]

                # match types
                # get ideal bond lengths, bondAB, bondCB
                # get ideal angle

Peter Eastman's avatar
Peter Eastman committed
2557
                if (type2 in types2 and ((type1 in types1 and type3 in types3) or (type3 in types1 and type1 in types3))): 
2558
2559
2560
                    if isConstrained:
                        hit = 1 
                    else:
Peter Eastman's avatar
Peter Eastman committed
2561
2562
2563
2564
                        hit = 1
                        bondAB = -1.0
                        bondCB = -1.0
                        swap = 0
2565
                        for bond in data.atomBonds[angle[1]]:
Peter Eastman's avatar
Peter Eastman committed
2566
2567
                            atom1 = data.bonds[bond].atom1
                            atom2 = data.bonds[bond].atom2
2568
                            length = data.bonds[bond].length
Peter Eastman's avatar
Peter Eastman committed
2569
                            if (atom1 == angle[0]):
2570
                                bondAB = length
Peter Eastman's avatar
Peter Eastman committed
2571
                            if (atom1 == angle[2]):
2572
                                bondCB = length
Peter Eastman's avatar
Peter Eastman committed
2573
                            if (atom2 == angle[2]):
2574
                                bondCB = length
Peter Eastman's avatar
Peter Eastman committed
2575
                            if (atom2 == angle[0]):
2576
2577
2578
2579
                                bondAB = length
                        
                        # check that ideal angle and bonds are set

Peter Eastman's avatar
Peter Eastman committed
2580
2581
2582
2583
                        if ('idealAngle' not in angleDict):
                           outputString = "AmoebaStretchBendGenerator: ideal angle is not set for following entry:\n"
                           outputString += "   %5d [%5d %5d %5d] [%5s %5s %5s]" % (count, angle[0], angle[1], angle[2], type1, type2, type3)
                           raise ValueError(outputString) 
2584

Peter Eastman's avatar
Peter Eastman committed
2585
2586
2587
2588
                        elif (bondAB < 0 or bondCB < 0):
                           outputString = "AmoebaStretchBendGenerator: bonds not set: %15.7e %15.7e. for following entry:" % (bondAB, bondCB)
                           outputString += "     %5d [%5d %5d %5d] [%5s %5s %5s]" % (count, angle[0], angle[1], angle[2], type1, type2, type3)
                           raise ValueError(outputString) 
2589
2590

                        else:
Peter Eastman's avatar
Peter Eastman committed
2591
                            if (verbose):
2592
                                 print "AmoebaStretchBendGenerator %5d [%5d %5d %5d] [%5s %5s %5s] %15.6f %15.6f %15.6f %15.6f" % (count, angle[0], angle[1], angle[2], type1, type2, type3, bondAB, bondCB, angleDict['idealAngle'], self.k1[i])
Peter Eastman's avatar
Peter Eastman committed
2593
                            force.addStretchBend(angle[0], angle[1], angle[2], bondAB, bondCB, angleDict['idealAngle']/radian, self.k1[i])
2594
2595
                    break

Peter Eastman's avatar
Peter Eastman committed
2596
            if (hit == 0 and verbose): 
2597
2598
2599
2600
2601
2602
2603
2604
2605
2606
2607
2608
2609
2610
                print "AmoebaStretchBendGenerator missing: %5d missing [%5d %5d %5d] [%5s %5s %5s] " % (count, angle[0], angle[1], angle[2], type1, type2, type3)

            count += 1

parsers["AmoebaStretchBendForce"] = AmoebaStretchBendGenerator.parseElement

#=============================================================================================

class AmoebaVdwGenerator:

    """A AmoebaVdwGenerator constructs a AmoebaVdwForce."""
    
    #=============================================================================================

Peter Eastman's avatar
Peter Eastman committed
2611
    def __init__(self, type, radiusrule, radiustype, radiussize, epsilonrule, vdw13Scale, vdw14Scale, vdw15Scale):
2612

Peter Eastman's avatar
Peter Eastman committed
2613
        self.type = type 
2614

Peter Eastman's avatar
Peter Eastman committed
2615
2616
2617
        self.radiusrule = radiusrule
        self.radiustype = radiustype
        self.radiussize = radiussize
2618

Peter Eastman's avatar
Peter Eastman committed
2619
        self.epsilonrule = epsilonrule
2620

Peter Eastman's avatar
Peter Eastman committed
2621
2622
2623
        self.vdw13Scale = vdw13Scale
        self.vdw14Scale = vdw14Scale
        self.vdw15Scale = vdw15Scale
2624

Peter Eastman's avatar
Peter Eastman committed
2625
        self.typeMap = {}
2626
2627
2628
2629
2630
2631
2632
2633
2634
2635

    #=============================================================================================

    @staticmethod
    def parseElement(element, forceField):

        # <AmoebaVdwForce type="BUFFERED-14-7" radiusrule="CUBIC-MEAN" radiustype="R-MIN" radiussize="DIAMETER" epsilonrule="HHG" vdw-13-scale="0.0" vdw-14-scale="1.0" vdw-15-scale="1.0" >
        #   <Vdw class="1" sigma="0.371" epsilon="0.46024" reduction="1.0" /> 
        #   <Vdw class="2" sigma="0.382" epsilon="0.422584" reduction="1.0" /> 
         
Peter Eastman's avatar
Peter Eastman committed
2636
2637
        generator = AmoebaVdwGenerator(element.attrib['type'], element.attrib['radiusrule'], element.attrib['radiustype'], element.attrib['radiussize'], element.attrib['epsilonrule'], 
                                        float(element.attrib['vdw-13-scale']), float(element.attrib['vdw-14-scale']), float(element.attrib['vdw-15-scale'])) 
2638
2639
2640
2641
2642
2643
2644
2645
2646
2647
        forceField._forces.append(generator)
        two_six = 1.122462048309372

        # types[] = [ sigma, epsilon, reductionFactor, class ]
        # sigma is modified based on radiustype and radiussize

        for atom in element.findall('Vdw'):
            types = forceField._findAtomTypes(atom, 1)
            if types is not None:

Peter Eastman's avatar
Peter Eastman committed
2648
                values = [float(atom.attrib['sigma']), float(atom.attrib['epsilon']), float(atom.attrib['reduction'])]
2649
2650
2651

                classType = atom.attrib['class']

Peter Eastman's avatar
Peter Eastman committed
2652
                if (generator.radiustype == 'SIGMA'):
2653
2654
                    values[0] *= two_six
      
Peter Eastman's avatar
Peter Eastman committed
2655
                if (generator.radiussize == 'DIAMETER'):
2656
2657
                    values[0] *= 0.5

Peter Eastman's avatar
Peter Eastman committed
2658
                values.append(classType)
2659
2660
2661
2662
2663

                for t in types[0]:
                    generator.typeMap[t] = values
    
            else:
Peter Eastman's avatar
Peter Eastman committed
2664
2665
                outputString = "AmoebaVdwGenerator: error getting type: %s" % (atom.attrib['class'])
                raise ValueError(outputString) 
2666
2667
2668
2669
2670
2671
2672
2673
    
    #=============================================================================================

    # Return a set containing the indices of particles bonded to particle with index=particleIndex

    #=============================================================================================

    @staticmethod
Peter Eastman's avatar
Peter Eastman committed
2674
    def getBondedParticleSet(particleIndex, data):
2675
2676
2677
2678
2679
2680

        bondedParticleSet = set()

        for bond in data.atomBonds[particleIndex]:
            atom1 = data.bonds[bond].atom1
            atom2 = data.bonds[bond].atom2
Peter Eastman's avatar
Peter Eastman committed
2681
2682
            if (atom1 != particleIndex):
                bondedParticleSet.add(atom1)
2683
            else:
Peter Eastman's avatar
Peter Eastman committed
2684
                bondedParticleSet.add(atom2)
2685
2686
2687
2688
2689

        return bondedParticleSet
          
    #=============================================================================================

Peter Eastman's avatar
Peter Eastman committed
2690
    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
2691

Peter Eastman's avatar
Peter Eastman committed
2692
        sigmaMap = {'ARITHMETIC':1, 'GEOMETRIC':1, 'CUBIC-MEAN':1}
2693
        epsilonMap = {'ARITHMETIC':1, 'GEOMETRIC':1, 'HARMONIC':1, 'HHG':1}
Peter Eastman's avatar
Peter Eastman committed
2694
        verbose = 0
2695
2696
2697
2698
2699
2700
2701
2702
2703
2704
2705

        # get or create force depending on whether it has already been added to the system

        existing = [sys.getForce(i) for i in range(sys.getNumForces())]
        existing = [f for f in existing if type(f) == mm.AmoebaVdwForce]
        if len(existing) == 0:
            force = mm.AmoebaVdwForce()
            sys.addForce(force)

            # sigma and epsilon combining rules

Peter Eastman's avatar
Peter Eastman committed
2706
            if ('sigmaCombiningRule' in args):
2707
                sigmaRule = args['sigmaCombiningRule'].upper()
Peter Eastman's avatar
Peter Eastman committed
2708
2709
                if (sigmaRule.upper() in sigmaMap):
                    force.setSigmaCombiningRule(sigmaRule.upper())
2710
                else:
Peter Eastman's avatar
Peter Eastman committed
2711
                    stringList = ' ' . join(str(x) for x in sigmaMap.keys())
2712
2713
                    print "sigma combining rule %s not recognized; valid values are %s; using default." % (sigmaRule, stringList)  
            else:
Peter Eastman's avatar
Peter Eastman committed
2714
                force.setSigmaCombiningRule(self.radiusrule)
2715

Peter Eastman's avatar
Peter Eastman committed
2716
            if ('epsilonCombiningRule' in args):
2717
                epsilonRule = args['epsilonCombiningRule'].upper()
Peter Eastman's avatar
Peter Eastman committed
2718
2719
                if (epsilonRule.upper() in epsilonMap):
                    force.setEpsilonCombiningRule(epsilonRule.upper())
2720
                else:
Peter Eastman's avatar
Peter Eastman committed
2721
                    stringList = ' ' . join(str(x) for x in epsilonMap.keys())
2722
2723
                    print "epsilon combining rule %s not recognized; valid values are %s; using default." % (epsilonRule, stringList)  
            else:
Peter Eastman's avatar
Peter Eastman committed
2724
                force.setEpsilonCombiningRule(self.epsilonrule)
2725
2726
2727
               
            # cutoff

Peter Eastman's avatar
Peter Eastman committed
2728
2729
            if ('vdwCutoff' in args):
                force.setCutoff(float(args['vdwCutoff']))
2730
            else:
Peter Eastman's avatar
Peter Eastman committed
2731
2732
2733
2734
                force.setCutoff(nonbondedCutoff)
            if (nonbondedMethod == PME):
                force.setPBC(1)
                force.setUseNeighborList(1)
2735
2736
2737
2738
2739
2740
2741
2742
2743
2744
2745
           
        else:
            force = existing[0]

        # add particles to force
        # throw error if particle type not available 

        for (i, atom) in enumerate(data.atoms):
            t = data.atomType[atom]
            if t in self.typeMap:

Peter Eastman's avatar
Peter Eastman committed
2746
2747
                values = self.typeMap[t]
                classIndex = int(values[3])
2748
2749
2750
   
                # ivIndex = index of bonded partner for hydrogens; otherwise ivIndex = particle index

Peter Eastman's avatar
Peter Eastman committed
2751
2752
2753
                ivIndex = i
                mass = sys.getParticleMass(i)/unit.dalton
                if (mass < 1.9 and len(data.atomBonds[i]) == 1):
2754
                    bondIndex = data.atomBonds[i][0]
Peter Eastman's avatar
Peter Eastman committed
2755
                    if (data.bonds[bondIndex].atom1 == i):
2756
2757
2758
2759
                        ivIndex = data.bonds[bondIndex].atom2
                    else:
                        ivIndex = data.bonds[bondIndex].atom1

Peter Eastman's avatar
Peter Eastman committed
2760
2761
                force.addParticle(ivIndex, classIndex, values[0], values[1], values[2])
                if (verbose):
2762
2763
2764
2765
2766
2767
2768
2769
2770
2771
2772
2773
                    print "Vdw %5d %5d %2d %15.7e %15.7e %15.7e" % (i, ivIndex, classIndex, values[0], values[1], values[2])
            else:
                raise ValueError('No vdw type for atom %s' % (atom.name))

        # set combining rules

        # set particle exclusions: self, 1-2 and 1-3 bonds
        # (1) collect in bondedParticleSets[i], 1-2 indices for all bonded partners of particle i
        # (2) add 1-2,1-3 and self to exclusion set

        bondedParticleSets = []
        for i in range(len(data.atoms)):
Peter Eastman's avatar
Peter Eastman committed
2774
            bondedParticleSets.append(AmoebaVdwGenerator.getBondedParticleSet(i, data))
2775
2776
2777
2778
2779
2780
2781
2782
2783

        for (i,atom) in enumerate(data.atoms):
 
            # 1-2 partners

            exclusionSet = bondedParticleSets[i].copy()

            # 1-3 partners

Peter Eastman's avatar
Peter Eastman committed
2784
            if (self.vdw13Scale == 0.0):
2785
                for bondedParticle in bondedParticleSets[i]:
Peter Eastman's avatar
Peter Eastman committed
2786
                    exclusionSet = exclusionSet.union(bondedParticleSets[bondedParticle])
2787
2788
2789
2790
2791

            # self

            exclusionSet.add(i)

Peter Eastman's avatar
Peter Eastman committed
2792
            if (verbose):
2793
                print "VdwExcl %5d %s || %s" % (i, str(exclusionSet), str(bondedParticleSets[i]))
Peter Eastman's avatar
Peter Eastman committed
2794
                for atomIndex in sorted(exclusionSet):
2795
                    atom = data.atoms[atomIndex]
Peter Eastman's avatar
Peter Eastman committed
2796
                    print "   %5d [%s %s %d]" % (atomIndex, atom.name, atom.residue.name, atom.residue.index)
2797
2798
                print "\n" 
            
Peter Eastman's avatar
Peter Eastman committed
2799
            force.setParticleExclusions(i, exclusionSet)
2800
2801
2802
2803
2804
2805
2806
2807
2808
2809
2810
2811
2812
2813
2814
2815
2816

parsers["AmoebaVdwForce"] = AmoebaVdwGenerator.parseElement

#=============================================================================================

class AmoebaMultipoleGenerator:

    #=============================================================================================

    """A AmoebaMultipoleGenerator constructs a AmoebaMultipoleForce."""
    
    #=============================================================================================

    def __init__(self, forceField,
                       direct11Scale, direct12Scale, direct13Scale, direct14Scale,
                       mpole12Scale,  mpole13Scale,  mpole14Scale,  mpole15Scale,
                       mutual11Scale, mutual12Scale, mutual13Scale, mutual14Scale,
Peter Eastman's avatar
Peter Eastman committed
2817
                       polar12Scale,  polar13Scale,  polar14Scale,  polar15Scale):
2818

Peter Eastman's avatar
Peter Eastman committed
2819
        self.forceField = forceField
2820

Peter Eastman's avatar
Peter Eastman committed
2821
2822
2823
2824
        self.direct11Scale = direct11Scale 
        self.direct12Scale = direct12Scale 
        self.direct13Scale = direct13Scale 
        self.direct14Scale = direct14Scale 
2825

Peter Eastman's avatar
Peter Eastman committed
2826
2827
2828
2829
        self.mpole12Scale = mpole12Scale 
        self.mpole13Scale = mpole13Scale 
        self.mpole14Scale = mpole14Scale 
        self.mpole15Scale = mpole15Scale 
2830

Peter Eastman's avatar
Peter Eastman committed
2831
2832
2833
2834
        self.mutual11Scale = mutual11Scale 
        self.mutual12Scale = mutual12Scale 
        self.mutual13Scale = mutual13Scale 
        self.mutual14Scale = mutual14Scale 
2835

Peter Eastman's avatar
Peter Eastman committed
2836
2837
2838
2839
        self.polar12Scale = polar12Scale 
        self.polar13Scale = polar13Scale 
        self.polar14Scale = polar14Scale 
        self.polar15Scale = polar15Scale 
2840

Peter Eastman's avatar
Peter Eastman committed
2841
2842
        self.typeMap = {}
        self.hasBeenCalled = 0
2843
2844
2845
2846
2847
2848

    #=============================================================================================
    # Set axis type
    #=============================================================================================

    @staticmethod
Peter Eastman's avatar
Peter Eastman committed
2849
    def setAxisType(kIndices):
2850
2851
2852

                # set axis type

Peter Eastman's avatar
Peter Eastman committed
2853
2854
2855
                kIndicesLen = len(kIndices)
                if (kIndicesLen > 3):
                    ky = kIndices[3]
2856
                else:
Peter Eastman's avatar
Peter Eastman committed
2857
                    ky = 0
2858
   
Peter Eastman's avatar
Peter Eastman committed
2859
2860
                if (kIndicesLen > 2):
                    kx = kIndices[2]
2861
                else:
Peter Eastman's avatar
Peter Eastman committed
2862
                    kx = 0
2863
   
Peter Eastman's avatar
Peter Eastman committed
2864
2865
                if (kIndicesLen > 1):
                    kz = kIndices[1]
2866
                else:
Peter Eastman's avatar
Peter Eastman committed
2867
                    kz = 0
2868

Peter Eastman's avatar
Peter Eastman committed
2869
2870
                while(len(kIndices) < 4):
                    kIndices.append(0)
2871
2872

                axisType = mm.AmoebaMultipoleForce.ZThenX
Peter Eastman's avatar
Peter Eastman committed
2873
                if (kz == 0):
2874
                    axisType = mm.AmoebaMultipoleForce.NoAxisType
Peter Eastman's avatar
Peter Eastman committed
2875
                if (kz != 0 and kx == 0):
2876
                    axisType = mm.AmoebaMultipoleForce.ZOnly
Peter Eastman's avatar
Peter Eastman committed
2877
                if (kz < 0 or kx < 0):
2878
                    axisType = mm.AmoebaMultipoleForce.Bisector
Peter Eastman's avatar
Peter Eastman committed
2879
                if (kx < 0 and ky < 0):
2880
                    axisType = mm.AmoebaMultipoleForce.ZBisect
Peter Eastman's avatar
Peter Eastman committed
2881
                if (kz < 0 and kx < 0 and ky  < 0):
2882
2883
                    axisType = mm.AmoebaMultipoleForce.ThreeFold

Peter Eastman's avatar
Peter Eastman committed
2884
2885
2886
                kIndices[1] = abs(kz)   
                kIndices[2] = abs(kx)   
                kIndices[3] = abs(ky)   
2887
2888
2889
2890
2891
2892
2893
2894
2895
2896
2897
2898

                return axisType

    #=============================================================================================

    @staticmethod
    def parseElement(element, forceField):

        #   <AmoebaMultipoleForce  direct11Scale="0.0"  direct12Scale="1.0"  direct13Scale="1.0"  direct14Scale="1.0"  mpole12Scale="0.0"  mpole13Scale="0.0"  mpole14Scale="0.4"  mpole15Scale="0.8"  mutual11Scale="1.0"  mutual12Scale="1.0"  mutual13Scale="1.0"  mutual14Scale="1.0"  polar12Scale="0.0"  polar13Scale="0.0"  polar14Intra="0.5"  polar14Scale="1.0"  polar15Scale="1.0"  > 
        # <Multipole class="1"    kz="2"    kx="4"    c0="-0.22620" d1="0.08214" d2="0.00000" d3="0.34883" q11="0.11775" q21="0.00000" q22="-1.02185" q31="-0.17555" q32="0.00000" q33="0.90410"  />
        # <Multipole class="2"    kz="1"    kx="3"    c0="-0.15245" d1="0.19517" d2="0.00000" d3="0.19687" q11="-0.20677" q21="0.00000" q22="-0.48084" q31="-0.01672" q32="0.00000" q33="0.68761"  />

Peter Eastman's avatar
Peter Eastman committed
2899
        generator = AmoebaMultipoleGenerator(forceField,
2900
2901
2902
2903
2904
2905
2906
2907
2908
2909
2910
2911
2912
2913
2914
2915
2916
2917
                                              element.attrib['direct11Scale'],
                                              element.attrib['direct12Scale'],
                                              element.attrib['direct13Scale'],
                                              element.attrib['direct14Scale'],

                                              element.attrib['mpole12Scale'],
                                              element.attrib['mpole13Scale'],
                                              element.attrib['mpole14Scale'],
                                              element.attrib['mpole15Scale'],

                                              element.attrib['mutual11Scale'],
                                              element.attrib['mutual12Scale'],
                                              element.attrib['mutual13Scale'],
                                              element.attrib['mutual14Scale'],

                                              element.attrib['polar12Scale'],
                                              element.attrib['polar13Scale'],
                                              element.attrib['polar14Scale'],
Peter Eastman's avatar
Peter Eastman committed
2918
                                              element.attrib['polar15Scale'])
2919
2920
2921
2922
2923
2924
2925
2926
2927
2928
2929



        forceField._forces.append(generator)

        # set type map: [ kIndices, multipoles, AMOEBA/OpenMM axis type]

        for atom in element.findall('Multipole'):
            types = forceField._findAtomTypes(atom, 1)
            if types is not None:

Peter Eastman's avatar
Peter Eastman committed
2930
                # print "Multipole Atom %s types=%s." %(atom.attrib['type'], str(types))
2931
2932
2933
2934
2935
2936
2937
2938

                # k-indices not provided default to 0

                kIndices = [int(atom.attrib['type'])]

                kStrings = [ 'kz', 'kx', 'ky' ]
                for kString in kStrings:
                    try:
Peter Eastman's avatar
Peter Eastman committed
2939
2940
                        if (atom.attrib[kString]):
                             kIndices.append(int(atom.attrib[kString]))
2941
2942
2943
2944
2945
                    except: 
                        pass

                # set axis type based on k-Indices 

Peter Eastman's avatar
Peter Eastman committed
2946
                axisType = AmoebaMultipoleGenerator.setAxisType(kIndices)
2947
2948
2949

                # set multipole

Peter Eastman's avatar
Peter Eastman committed
2950
                charge = float(atom.attrib['c0'])
2951
 
Peter Eastman's avatar
Peter Eastman committed
2952
2953
2954
2955
2956
2957
2958
2959
2960
2961
2962
2963
2964
                conversion = 1.0
                dipole = [ conversion*float(atom.attrib['d1']), conversion*float(atom.attrib['d2']), conversion*float(atom.attrib['d3'])]

                quadrupole = []
                quadrupole.append(conversion*float(atom.attrib['q11']))
                quadrupole.append(conversion*float(atom.attrib['q21']))
                quadrupole.append(conversion*float(atom.attrib['q31']))
                quadrupole.append(conversion*float(atom.attrib['q21']))
                quadrupole.append(conversion*float(atom.attrib['q22']))
                quadrupole.append(conversion*float(atom.attrib['q32']))
                quadrupole.append(conversion*float(atom.attrib['q31']))
                quadrupole.append(conversion*float(atom.attrib['q32']))
                quadrupole.append(conversion*float(atom.attrib['q33']))
2965
2966

                for t in types[0]:
Peter Eastman's avatar
Peter Eastman committed
2967
                    if (t not in generator.typeMap):
2968
2969
                        generator.typeMap[t] = []

Peter Eastman's avatar
Peter Eastman committed
2970
2971
2972
2973
2974
2975
2976
2977
                    valueMap = dict()
                    valueMap['classIndex'] = atom.attrib['type']
                    valueMap['kIndices'] = kIndices
                    valueMap['charge'] = charge 
                    valueMap['dipole'] = dipole
                    valueMap['quadrupole'] = quadrupole
                    valueMap['axisType'] = axisType
                    generator.typeMap[t].append(valueMap)
2978
2979
    
            else:
Peter Eastman's avatar
Peter Eastman committed
2980
2981
                outputString = "AmoebaMultipoleGenerator: error getting type for multipole: %s" % (atom.attrib['class'])
                raise ValueError(outputString) 
2982
2983
2984
2985
2986
2987
2988
    
        # polarization parameters
 
        for atom in element.findall('Polarize'):
            types = forceField._findAtomTypes(atom, 1)
            if types is not None:

Peter Eastman's avatar
Peter Eastman committed
2989
2990
2991
2992
                classIndex = atom.attrib['type']
                polarizability = float(atom.attrib['polarizability'])
                thole = float(atom.attrib['thole'])
                if (thole == 0):
2993
2994
                    pdamp = 0
                else:
Peter Eastman's avatar
Peter Eastman committed
2995
                    pdamp = pow(polarizability, 1.0/6.0)
2996

Peter Eastman's avatar
Peter Eastman committed
2997
2998
                pgrpMap = dict()
                for index in range(1, 7):
2999
                    pgrp = 'pgrp' + str(index)
Peter Eastman's avatar
Peter Eastman committed
3000
                    if (pgrp in atom.attrib):
3001
3002
3003
                        pgrpMap[int(atom.attrib[pgrp])] = -1

                for t in types[0]:
Peter Eastman's avatar
Peter Eastman committed
3004
3005
3006
                    if (t not in generator.typeMap):
                        outputString = "AmoebaMultipoleGenerator: polarize type not present: %s" % (atom.attrib['type'])
                        raise ValueError(outputString) 
3007
3008
                    else:
                        typeMapList = generator.typeMap[t]
Peter Eastman's avatar
Peter Eastman committed
3009
3010
3011
3012
3013
3014
3015
3016
3017
3018
3019
3020
3021
3022
3023
                        hit = 0
                        for (ii, typeMap) in enumerate(typeMapList):

                            if (typeMap['classIndex'] == classIndex):
                                typeMap['polarizability'] = polarizability
                                typeMap['thole'] = thole
                                typeMap['pdamp'] = pdamp
                                typeMap['pgrpMap'] = pgrpMap 
                                typeMapList[ii] = typeMap
                                hit = 1
                                #print "Adding polarize for %s map=%d len=%d keys=%s" % (classIndex, ii, len(typeMapList), str(typeMap.keys()))

                        if (hit == 0):
                            outputString = "AmoebaMultipoleGenerator: error getting type for polarize: class index=%s not in multipole list?" % (atom.attrib['class'])
                            raise ValueError(outputString) 
3024
3025
    
            else:
Peter Eastman's avatar
Peter Eastman committed
3026
3027
                outputString = "AmoebaMultipoleGenerator: error getting type for polarize: %s" % (atom.attrib['class'])
                raise ValueError(outputString) 
3028
3029
3030
    
    #=============================================================================================

Peter Eastman's avatar
Peter Eastman committed
3031
    def setPolarGroups(self, data, bonded12ParticleSets, force):
3032
3033
3034
3035
3036

        for (atomIndex, atom) in enumerate(data.atoms):

            # assign multipole parameters via only 1-2 connected atoms

Peter Eastman's avatar
Peter Eastman committed
3037
3038
3039
3040
            multipoleDict = atom.multipoleDict
            pgrpMap = multipoleDict['pgrpMap']
            bondedAtomIndices = bonded12ParticleSets[atomIndex]
            atom.stage = -1
3041
3042
3043
3044
            atom.polarizationGroupSet = list()
            atom.polarizationGroups[atomIndex] = 1
            for bondedAtomIndex in bondedAtomIndices:
                bondedAtomType = int(data.atomType[data.atoms[bondedAtomIndex]])
Peter Eastman's avatar
Peter Eastman committed
3045
3046
                bondedAtom = data.atoms[bondedAtomIndex]
                if (bondedAtomType in pgrpMap):
3047
3048
3049
3050
3051
3052
3053
                    atom.polarizationGroups[bondedAtomIndex] = 1
                    bondedAtom.polarizationGroups[atomIndex] = 1
                    
        # pgrp11

        for (atomIndex, atom) in enumerate(data.atoms):

Peter Eastman's avatar
Peter Eastman committed
3054
            if (len( data.atoms[atomIndex].polarizationGroupSet) > 0):
3055
3056
                continue

Peter Eastman's avatar
Peter Eastman committed
3057
3058
            group = set()
            visited = set()
3059
3060
            notVisited = set()
            for pgrpAtomIndex in atom.polarizationGroups:
Peter Eastman's avatar
Peter Eastman committed
3061
3062
3063
3064
                group.add(pgrpAtomIndex)
                notVisited.add(pgrpAtomIndex)
            visited.add(atomIndex)
            while(len(notVisited) > 0):
3065
                nextAtom = notVisited.pop()
Peter Eastman's avatar
Peter Eastman committed
3066
3067
                if (nextAtom not in visited):
                   visited.add(nextAtom)
3068
                   for ii in data.atoms[nextAtom].polarizationGroups:
Peter Eastman's avatar
Peter Eastman committed
3069
3070
3071
                       group.add(ii)
                       if (ii not in visited):
                           notVisited.add(ii)
3072
3073
3074

            pGroup = group
            for pgrpAtomIndex in group:
Peter Eastman's avatar
Peter Eastman committed
3075
                data.atoms[pgrpAtomIndex].polarizationGroupSet.append(pGroup)
3076
3077

        for (atomIndex, atom) in enumerate(data.atoms):
Peter Eastman's avatar
Peter Eastman committed
3078
3079
            atom.polarizationGroupSet[0] = sorted(atom.polarizationGroupSet[0])
            force.setCovalentMap(atomIndex, mm.AmoebaMultipoleForce.PolarizationCovalent11, atom.polarizationGroupSet[0])
3080
3081
3082
3083
3084

        # pgrp12

        for (atomIndex, atom) in enumerate(data.atoms):

Peter Eastman's avatar
Peter Eastman committed
3085
            if (len( data.atoms[atomIndex].polarizationGroupSet) > 1):
3086
3087
                continue

Peter Eastman's avatar
Peter Eastman committed
3088
            pgrp11 = set(atom.polarizationGroupSet[0])
3089
3090
3091
            pgrp12 = set()
            for pgrpAtomIndex in pgrp11:
                for bonded12 in bonded12ParticleSets[pgrpAtomIndex]:
Peter Eastman's avatar
Peter Eastman committed
3092
                    pgrp12 = pgrp12.union(data.atoms[bonded12].polarizationGroupSet[0])
3093
3094
            pgrp12 = pgrp12 - pgrp11
            for pgrpAtomIndex in pgrp11:
Peter Eastman's avatar
Peter Eastman committed
3095
                data.atoms[pgrpAtomIndex].polarizationGroupSet.append(pgrp12)
3096
3097
                
        for (atomIndex, atom) in enumerate(data.atoms):
Peter Eastman's avatar
Peter Eastman committed
3098
3099
            atom.polarizationGroupSet[1] = sorted(atom.polarizationGroupSet[1])
            force.setCovalentMap(atomIndex, mm.AmoebaMultipoleForce.PolarizationCovalent12, atom.polarizationGroupSet[1])
3100
3101
3102
3103
3104

        # pgrp13

        for (atomIndex, atom) in enumerate(data.atoms):

Peter Eastman's avatar
Peter Eastman committed
3105
            if (len(data.atoms[atomIndex].polarizationGroupSet) > 2):
3106
3107
                continue

Peter Eastman's avatar
Peter Eastman committed
3108
3109
            pgrp11 = set(atom.polarizationGroupSet[0])
            pgrp12 = set(atom.polarizationGroupSet[1])
3110
3111
3112
            pgrp13 = set()
            for pgrpAtomIndex in pgrp12:
                for bonded12 in bonded12ParticleSets[pgrpAtomIndex]:
Peter Eastman's avatar
Peter Eastman committed
3113
                    pgrp13 = pgrp13.union(data.atoms[bonded12].polarizationGroupSet[0])
3114
            pgrp13 = pgrp13 - pgrp12
Peter Eastman's avatar
Peter Eastman committed
3115
            pgrp13 = pgrp13 - set(pgrp11)
3116
            for pgrpAtomIndex in pgrp11:
Peter Eastman's avatar
Peter Eastman committed
3117
                data.atoms[pgrpAtomIndex].polarizationGroupSet.append(pgrp13)
3118
3119
                
        for (atomIndex, atom) in enumerate(data.atoms):
Peter Eastman's avatar
Peter Eastman committed
3120
3121
            atom.polarizationGroupSet[2] = sorted(atom.polarizationGroupSet[2])
            force.setCovalentMap(atomIndex, mm.AmoebaMultipoleForce.PolarizationCovalent13, atom.polarizationGroupSet[2])
3122
3123
3124
3125
3126

        # pgrp14

        for (atomIndex, atom) in enumerate(data.atoms):

Peter Eastman's avatar
Peter Eastman committed
3127
            if (len(data.atoms[atomIndex].polarizationGroupSet) > 3):
3128
3129
                continue

Peter Eastman's avatar
Peter Eastman committed
3130
3131
3132
            pgrp11 = set(atom.polarizationGroupSet[0])
            pgrp12 = set(atom.polarizationGroupSet[1])
            pgrp13 = set(atom.polarizationGroupSet[2])
3133
3134
3135
            pgrp14 = set()
            for pgrpAtomIndex in pgrp13:
                for bonded12 in bonded12ParticleSets[pgrpAtomIndex]:
Peter Eastman's avatar
Peter Eastman committed
3136
                    pgrp14 = pgrp14.union(data.atoms[bonded12].polarizationGroupSet[0])
3137
3138
3139

            pgrp14 = pgrp14 - pgrp13
            pgrp14 = pgrp14 - pgrp12
Peter Eastman's avatar
Peter Eastman committed
3140
            pgrp14 = pgrp14 - set(pgrp11)
3141
3142

            for pgrpAtomIndex in pgrp11:
Peter Eastman's avatar
Peter Eastman committed
3143
                data.atoms[pgrpAtomIndex].polarizationGroupSet.append(pgrp14)
3144
3145
                
        for (atomIndex, atom) in enumerate(data.atoms):
Peter Eastman's avatar
Peter Eastman committed
3146
3147
            atom.polarizationGroupSet[3] = sorted(atom.polarizationGroupSet[3])
            force.setCovalentMap(atomIndex, mm.AmoebaMultipoleForce.PolarizationCovalent14, atom.polarizationGroupSet[3])
3148
3149
3150

    #=============================================================================================

Peter Eastman's avatar
Peter Eastman committed
3151
    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
3152

Peter Eastman's avatar
Peter Eastman committed
3153
        if (self.hasBeenCalled ):
3154
3155
3156
3157
3158
3159
3160
3161
3162
3163
3164
3165
3166
3167
            return
        self.hasBeenCalled = 1

        methodMap = {NoCutoff:mm.AmoebaMultipoleForce.NoCutoff,
                     PME:mm.AmoebaMultipoleForce.PME}
        verbose = 0

        # get or create force depending on whether it has already been added to the system

        existing = [sys.getForce(i) for i in range(sys.getNumForces())]
        existing = [f for f in existing if type(f) == mm.AmoebaMultipoleForce]
        if len(existing) == 0:
            force = mm.AmoebaMultipoleForce()
            sys.addForce(force)
Peter Eastman's avatar
Peter Eastman committed
3168
            if (nonbondedMethod not in methodMap): 
3169
                print "Warning: AmoebaMultipoleForce: cutoff method not available using NoCutoff."
Peter Eastman's avatar
Peter Eastman committed
3170
                force.setNonbondedMethod(mm.AmoebaMultipoleForce.NoCutoff)
3171
            else:
Peter Eastman's avatar
Peter Eastman committed
3172
                force.setNonbondedMethod(methodMap[nonbondedMethod])
3173
3174
            force.setCutoffDistance(nonbondedCutoff)

Peter Eastman's avatar
Peter Eastman committed
3175
            if ('ewaldErrorTolerance' in args):
3176
3177
                force.setEwaldErrorTolerance(float(args['ewaldErrorTolerance']))

Peter Eastman's avatar
Peter Eastman committed
3178
            if ('polarization' in args):
3179
                polarizationType = args['polarization']
Peter Eastman's avatar
Peter Eastman committed
3180
                if (polarizationType.lower() == 'direct'):
3181
3182
3183
3184
                    force.setPolarizationType(mm.AmoebaMultipoleForce.Direct)
                else:
                    force.setPolarizationType(mm.AmoebaMultipoleForce.Mutual)

Peter Eastman's avatar
Peter Eastman committed
3185
3186
            if ('aEwald' in args):
                force.setAEwald(float(args['aEwald']))
3187

Peter Eastman's avatar
Peter Eastman committed
3188
            if ('pmeGridDimensions' in args):
3189
3190
                force.setPmeGridDimensions(args['pmeGridDimensions'])

Peter Eastman's avatar
Peter Eastman committed
3191
3192
            if ('mutualInducedMaxIterations' in args):
                force.setMutualInducedMaxIterations(int(args['mutualInducedMaxIterations']))
3193

Peter Eastman's avatar
Peter Eastman committed
3194
            if ('mutualInducedTargetEpsilon' in args):
3195
3196
3197
3198
3199
3200
3201
3202
3203
3204
3205
3206
3207
3208
                force.setMutualInducedTargetEpsilon(float(args['mutualInducedTargetEpsilon']))

        else:
            force = existing[0]

        # add particles to force
        # throw error if particle type not available 

        # get 1-2, 1-3, 1-4, 1-5 bonded sets

        # 1-2

        bonded12ParticleSets = []
        for i in range(len(data.atoms)):
Peter Eastman's avatar
Peter Eastman committed
3209
3210
3211
            bonded12ParticleSet = AmoebaVdwGenerator.getBondedParticleSet(i, data)
            bonded12ParticleSet = set(sorted(bonded12ParticleSet))
            bonded12ParticleSets.append(bonded12ParticleSet)
3212
3213
3214
3215
3216

        # 1-3

        bonded13ParticleSets = []
        for i in range(len(data.atoms)):
Peter Eastman's avatar
Peter Eastman committed
3217
            bonded13Set = set()
3218
3219
            bonded12ParticleSet = bonded12ParticleSets[i]
            for j in bonded12ParticleSet: 
Peter Eastman's avatar
Peter Eastman committed
3220
                bonded13Set = bonded13Set.union(bonded12ParticleSets[j])
3221
3222
3223
3224

            # remove 1-2 and self from set

            bonded13Set = bonded13Set - bonded12ParticleSet
Peter Eastman's avatar
Peter Eastman committed
3225
            selfSet = set()
3226
3227
            selfSet.add(i)
            bonded13Set = bonded13Set - selfSet
Peter Eastman's avatar
Peter Eastman committed
3228
3229
            bonded13Set = set(sorted(bonded13Set))
            bonded13ParticleSets.append(bonded13Set)
3230
3231
3232
3233
3234

        # 1-4

        bonded14ParticleSets = []
        for i in range(len(data.atoms)):
Peter Eastman's avatar
Peter Eastman committed
3235
3236
            bonded14Set = set()
            bonded13ParticleSet = bonded13ParticleSets[i]
3237
            for j in bonded13ParticleSet: 
Peter Eastman's avatar
Peter Eastman committed
3238
                bonded14Set = bonded14Set.union(bonded12ParticleSets[j])
3239
3240
3241
3242
3243
           
            # remove 1-3, 1-2 and self from set

            bonded14Set = bonded14Set - bonded12ParticleSets[i]
            bonded14Set = bonded14Set - bonded13ParticleSet
Peter Eastman's avatar
Peter Eastman committed
3244
            selfSet = set()
3245
3246
            selfSet.add(i)
            bonded14Set = bonded14Set - selfSet
Peter Eastman's avatar
Peter Eastman committed
3247
3248
            bonded14Set = set(sorted(bonded14Set))
            bonded14ParticleSets.append(bonded14Set)
3249
3250
3251
3252
3253

        # 1-5

        bonded15ParticleSets = []
        for i in range(len(data.atoms)):
Peter Eastman's avatar
Peter Eastman committed
3254
3255
            bonded15Set = set()
            bonded14ParticleSet = bonded14ParticleSets[i]
3256
            for j in bonded14ParticleSet: 
Peter Eastman's avatar
Peter Eastman committed
3257
                bonded15Set = bonded15Set.union(bonded12ParticleSets[j])
3258
3259
3260
3261
3262
3263

            # remove 1-4, 1-3, 1-2 and self from set

            bonded15Set = bonded15Set - bonded12ParticleSets[i]
            bonded15Set = bonded15Set - bonded13ParticleSets[i]
            bonded15Set = bonded15Set - bonded14ParticleSet
Peter Eastman's avatar
Peter Eastman committed
3264
            selfSet = set()
3265
3266
            selfSet.add(i)
            bonded15Set = bonded15Set - selfSet
Peter Eastman's avatar
Peter Eastman committed
3267
3268
            bonded15Set = set(sorted(bonded15Set))
            bonded15ParticleSets.append(bonded15Set)
3269
3270
3271
3272
3273

        for (atomIndex, atom) in enumerate(data.atoms):
            t = data.atomType[atom]
            if t in self.typeMap:

Peter Eastman's avatar
Peter Eastman committed
3274
3275
                multipoleList = self.typeMap[t]
                hit = 0
3276
3277
3278
3279
3280
3281
                savedMultipoleDict = 0

                # assign multipole parameters via only 1-2 connected atoms

                for multipoleDict in multipoleList:

Peter Eastman's avatar
Peter Eastman committed
3282
                    if (hit != 0):
3283
3284
                        break

Peter Eastman's avatar
Peter Eastman committed
3285
                    kIndices = multipoleDict['kIndices']
3286
    
Peter Eastman's avatar
Peter Eastman committed
3287
3288
3289
                    kz = kIndices[1]   
                    kx = kIndices[2]
                    ky = kIndices[3]
3290
3291
3292
3293
3294
3295

                    # assign multipole parameters
                    #    (1) get bonded partners
                    #    (2) match parameter types
    
                    bondedAtomIndices = bonded12ParticleSets[atomIndex]
Peter Eastman's avatar
Peter Eastman committed
3296
3297
3298
                    zaxis = -1
                    xaxis = -1
                    yaxis = -1
3299
3300
                    for bondedAtomZIndex in bondedAtomIndices:

Peter Eastman's avatar
Peter Eastman committed
3301
                       if (hit != 0):
3302
3303
3304
                           break

                       bondedAtomZType = int(data.atomType[data.atoms[bondedAtomZIndex]])
Peter Eastman's avatar
Peter Eastman committed
3305
3306
                       bondedAtomZ = data.atoms[bondedAtomZIndex]
                       if (bondedAtomZType == kz):
3307
                          for bondedAtomXIndex in bondedAtomIndices:
Peter Eastman's avatar
Peter Eastman committed
3308
                              if (bondedAtomXIndex == bondedAtomZIndex or hit != 0):
3309
3310
                                  continue
                              bondedAtomXType = int(data.atomType[data.atoms[bondedAtomXIndex]])
Peter Eastman's avatar
Peter Eastman committed
3311
3312
3313
3314
                              if (bondedAtomXType == kx):
                                  if (ky == 0):
                                      zaxis = bondedAtomZIndex
                                      xaxis = bondedAtomXIndex
3315
                                      savedMultipoleDict = multipoleDict
Peter Eastman's avatar
Peter Eastman committed
3316
                                      hit = 1
3317
3318
                                  else:
                                      for bondedAtomYIndex in bondedAtomIndices:
Peter Eastman's avatar
Peter Eastman committed
3319
                                          if (bondedAtomYIndex == bondedAtomZIndex or bondedAtomYIndex == bondedAtomXIndex or hit != 0):
3320
3321
                                              continue
                                          bondedAtomYType = int(data.atomType[data.atoms[bondedAtomYIndex]])
Peter Eastman's avatar
Peter Eastman committed
3322
3323
3324
3325
                                          if (bondedAtomYType == ky):
                                              zaxis = bondedAtomZIndex
                                              xaxis = bondedAtomXIndex
                                              yaxis = bondedAtomYIndex
3326
                                              savedMultipoleDict = multipoleDict
Peter Eastman's avatar
Peter Eastman committed
3327
                                              hit = 2
3328
3329
3330
3331
3332
                                         
                # assign multipole parameters via 1-2 and 1-3 connected atoms

                for multipoleDict in multipoleList:

Peter Eastman's avatar
Peter Eastman committed
3333
                    if (hit != 0):
3334
3335
                        break

Peter Eastman's avatar
Peter Eastman committed
3336
                    kIndices = multipoleDict['kIndices']
3337
    
Peter Eastman's avatar
Peter Eastman committed
3338
3339
3340
                    kz = kIndices[1]   
                    kx = kIndices[2]
                    ky = kIndices[3]
3341
3342
3343
3344
3345
3346
3347
3348
    
                    # assign multipole parameters
                    #    (1) get bonded partners
                    #    (2) match parameter types
    
                    bondedAtom12Indices = bonded12ParticleSets[atomIndex]
                    bondedAtom13Indices = bonded13ParticleSets[atomIndex]

Peter Eastman's avatar
Peter Eastman committed
3349
3350
3351
                    zaxis = -1
                    xaxis = -1
                    yaxis = -1
3352
3353
3354

                    for bondedAtomZIndex in bondedAtom12Indices:

Peter Eastman's avatar
Peter Eastman committed
3355
                       if (hit != 0):
3356
3357
3358
                           break

                       bondedAtomZType = int(data.atomType[data.atoms[bondedAtomZIndex]])
Peter Eastman's avatar
Peter Eastman committed
3359
                       bondedAtomZ = data.atoms[bondedAtomZIndex]
3360

Peter Eastman's avatar
Peter Eastman committed
3361
                       if (bondedAtomZType == kz):
3362
3363
                          for bondedAtomXIndex in bondedAtom13Indices:

Peter Eastman's avatar
Peter Eastman committed
3364
                              if (bondedAtomXIndex == bondedAtomZIndex or hit != 0):
3365
3366
                                  continue
                              bondedAtomXType = int(data.atomType[data.atoms[bondedAtomXIndex]])
Peter Eastman's avatar
Peter Eastman committed
3367
3368
3369
3370
                              if (bondedAtomXType == kx and bondedAtomZIndex in bonded12ParticleSets[bondedAtomXIndex]):
                                  if (ky == 0):
                                      zaxis = bondedAtomZIndex
                                      xaxis = bondedAtomXIndex
3371
                                      savedMultipoleDict = multipoleDict
Peter Eastman's avatar
Peter Eastman committed
3372
                                      hit = 3
3373
3374
                                  else:
                                      for bondedAtomYIndex in bondedAtom13Indices:
Peter Eastman's avatar
Peter Eastman committed
3375
                                          if (bondedAtomYIndex == bondedAtomZIndex or bondedAtomYIndex == bondedAtomXIndex or hit != 0):
3376
3377
                                              continue
                                          bondedAtomYType = int(data.atomType[data.atoms[bondedAtomYIndex]])
Peter Eastman's avatar
Peter Eastman committed
3378
3379
3380
3381
                                          if (bondedAtomYType == ky and bondedAtomZIndex in bonded12ParticleSets[bondedAtomYIndex]):
                                              zaxis = bondedAtomZIndex
                                              xaxis = bondedAtomXIndex
                                              yaxis = bondedAtomYIndex
3382
                                              savedMultipoleDict = multipoleDict
Peter Eastman's avatar
Peter Eastman committed
3383
                                              hit = 4
3384
3385
3386
3387
3388
                                         
                # assign multipole parameters via only a z-defining atom

                for multipoleDict in multipoleList:

Peter Eastman's avatar
Peter Eastman committed
3389
                    if (hit != 0):
3390
3391
                        break

Peter Eastman's avatar
Peter Eastman committed
3392
                    kIndices = multipoleDict['kIndices']
3393
    
Peter Eastman's avatar
Peter Eastman committed
3394
3395
                    kz = kIndices[1]   
                    kx = kIndices[2]   
3396
    
Peter Eastman's avatar
Peter Eastman committed
3397
3398
3399
                    zaxis = -1
                    xaxis = -1
                    yaxis = -1
3400
3401
3402

                    for bondedAtomZIndex in bondedAtom12Indices:

Peter Eastman's avatar
Peter Eastman committed
3403
                        if (hit != 0):
3404
3405
3406
                            break

                        bondedAtomZType = int(data.atomType[data.atoms[bondedAtomZIndex]])
Peter Eastman's avatar
Peter Eastman committed
3407
                        bondedAtomZ = data.atoms[bondedAtomZIndex]
3408

Peter Eastman's avatar
Peter Eastman committed
3409
3410
                        if (kx == 0 and kz == bondedAtomZType):
                            kz = bondedAtomZIndex
3411
                            savedMultipoleDict = multipoleDict
Peter Eastman's avatar
Peter Eastman committed
3412
                            hit = 5
3413
3414
3415
3416
3417

                # assign multipole parameters via no connected atoms

                for multipoleDict in multipoleList:

Peter Eastman's avatar
Peter Eastman committed
3418
                    if (hit != 0):
3419
3420
                        break

Peter Eastman's avatar
Peter Eastman committed
3421
                    kIndices = multipoleDict['kIndices']
3422
    
Peter Eastman's avatar
Peter Eastman committed
3423
                    kz = kIndices[1]   
3424
    
Peter Eastman's avatar
Peter Eastman committed
3425
3426
3427
                    zaxis = -1
                    xaxis = -1
                    yaxis = -1
3428

Peter Eastman's avatar
Peter Eastman committed
3429
                    if (kz == 0):
3430
                        savedMultipoleDict = multipoleDict
Peter Eastman's avatar
Peter Eastman committed
3431
                        hit = 6
3432
3433
3434
                
                # add particle if there was a hit

Peter Eastman's avatar
Peter Eastman committed
3435
3436
3437
3438
                if (hit != 0):
                    if (verbose):
                        print "Multipole hit for Atom %5d %4s of %4s %4d type=%5s hitBranch=%d axes: %4d %4d %4d axisType=%1d." % (atomIndex,
                               atom.name, atom.residue.name, atom.residue.index, t, hit, zaxis, xaxis, yaxis, savedMultipoleDict['axisType'])
3439

Peter Eastman's avatar
Peter Eastman committed
3440
                    atom.multipoleDict = savedMultipoleDict
3441
                    atom.polarizationGroups = dict()
Peter Eastman's avatar
Peter Eastman committed
3442
                    newIndex = force.addParticle(savedMultipoleDict['charge'], savedMultipoleDict['dipole'], savedMultipoleDict['quadrupole'], savedMultipoleDict['axisType'],
3443
                                                                 zaxis, xaxis, yaxis, savedMultipoleDict['thole'], savedMultipoleDict['pdamp'], savedMultipoleDict['polarizability'])
Peter Eastman's avatar
Peter Eastman committed
3444
3445
3446
3447
3448
                    if (atomIndex == newIndex):
                        force.setCovalentMap(atomIndex, mm.AmoebaMultipoleForce.Covalent12, bonded12ParticleSets[atomIndex])
                        force.setCovalentMap(atomIndex, mm.AmoebaMultipoleForce.Covalent13, bonded13ParticleSets[atomIndex])
                        force.setCovalentMap(atomIndex, mm.AmoebaMultipoleForce.Covalent14, bonded14ParticleSets[atomIndex])
                        force.setCovalentMap(atomIndex, mm.AmoebaMultipoleForce.Covalent15, bonded15ParticleSets[atomIndex])
3449
                    else:
Peter Eastman's avatar
Peter Eastman committed
3450
                        raise ValueError("Atom %s of %s %d is out of synch!." %(atom.name, atom.residue.name, atom.residue.index))
3451
                else:
Peter Eastman's avatar
Peter Eastman committed
3452
                    raise ValueError("Atom %s of %s %d was not assigned." %(atom.name, atom.residue.name, atom.residue.index))
3453
            else:
Peter Eastman's avatar
Peter Eastman committed
3454
                raise ValueError('No multipole type for atom %s %s %d' % (atom.name, atom.residue.name, atom.residue.index))
3455
3456
3457

        # set polar groups

Peter Eastman's avatar
Peter Eastman committed
3458
        self.setPolarGroups(data, bonded12ParticleSets, force)
3459
3460
3461
3462
3463
3464
3465
3466
3467
3468
3469

parsers["AmoebaMultipoleForce"] = AmoebaMultipoleGenerator.parseElement

#=============================================================================================

class AmoebaWcaDispersionGenerator:

    """A AmoebaWcaDispersionGenerator constructs a AmoebaWcaDispersionForce."""
    
    #=========================================================================================

Peter Eastman's avatar
Peter Eastman committed
3470
    def __init__(self, epso, epsh, rmino, rminh, awater, slevy, dispoff, shctd):
3471

Peter Eastman's avatar
Peter Eastman committed
3472
3473
3474
3475
3476
3477
3478
3479
        self.epso = epso 
        self.epsh = epsh 
        self.rmino = rmino
        self.rminh = rminh
        self.awater = awater
        self.slevy = slevy
        self.dispoff = dispoff
        self.shctd = shctd 
3480

Peter Eastman's avatar
Peter Eastman committed
3481
        self.typeMap = {}
3482
3483
3484
3485
3486
3487
3488
3489
3490
3491

    #=========================================================================================

    @staticmethod
    def parseElement(element, forceField):

        #  <AmoebaWcaDispersionForce epso="0.46024" epsh="0.056484" rmino="0.17025" rminh="0.13275" awater="33.428" slevy="1.0"  dispoff="0.026" shctd="0.81" >
        #   <WcaDispersion class="1" radius="0.1855" epsilon="0.46024" />
        #   <WcaDispersion class="2" radius="0.191" epsilon="0.422584" />
      
Peter Eastman's avatar
Peter Eastman committed
3492
        generator = AmoebaWcaDispersionGenerator(element.attrib['epso'],
3493
3494
3495
3496
3497
3498
                                                  element.attrib['epsh'],
                                                  element.attrib['rmino'],
                                                  element.attrib['rminh'],
                                                  element.attrib['awater'], 
                                                  element.attrib['slevy'],
                                                  element.attrib['dispoff'],
Peter Eastman's avatar
Peter Eastman committed
3499
                                                  element.attrib['shctd']) 
3500
3501
3502
3503
3504
3505
3506
3507
3508
3509
3510
3511
        forceField._forces.append(generator)

        # typeMap[] = [ radius, epsilon ]

        for atom in element.findall('WcaDispersion'):
            types = forceField._findAtomTypes(atom, 1)
            if types is not None:

                values = [float(atom.attrib['radius']), float(atom.attrib['epsilon'])]
                for t in types[0]:
                    generator.typeMap[t] = values
            else:
Peter Eastman's avatar
Peter Eastman committed
3512
3513
                outputString = "AmoebaWcaDispersionGenerator: error getting type: %s" % (atom.attrib['class'])
                raise ValueError(outputString) 
3514
3515
3516
    
    #=========================================================================================
    
Peter Eastman's avatar
Peter Eastman committed
3517
    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
3518
3519
3520
3521
3522
3523
3524
3525
3526
3527
3528
3529
3530
3531
3532
3533

        verbose = 0

        # get or create force depending on whether it has already been added to the system

        existing = [sys.getForce(i) for i in range(sys.getNumForces())]
        existing = [f for f in existing if type(f) == mm.AmoebaWcaDispersionForce]
        if len(existing) == 0:
            force = mm.AmoebaWcaDispersionForce()
            sys.addForce(force)
        else:
            force = existing[0]

        # add particles to force
        # throw error if particle type not available 

Peter Eastman's avatar
Peter Eastman committed
3534
3535
3536
3537
3538
3539
3540
3541
        force.setEpso(   float(self.epso   ))
        force.setEpsh(   float(self.epsh   ))
        force.setRmino(  float(self.rmino  ))
        force.setRminh(  float(self.rminh  ))
        force.setDispoff(float(self.dispoff))
        force.setSlevy(  float(self.slevy  ))
        force.setAwater( float(self.awater ))
        force.setShctd(  float(self.shctd  ))
3542
3543
3544
3545
3546

        for (i, atom) in enumerate(data.atoms):
            t = data.atomType[atom]
            if t in self.typeMap:

Peter Eastman's avatar
Peter Eastman committed
3547
3548
3549
                values = self.typeMap[t]
                force.addParticle(values[0], values[1])
                if (verbose):
3550
3551
3552
3553
3554
3555
3556
3557
3558
3559
3560
3561
3562
3563
                    print "WcaDispersion %5d %15.7e %15.7e" % (i, values[0], values[1])
            else:
                raise ValueError('No WcaDispersion type for atom %s of %s %d' % (atom.name, atom.residue.name, atom.residue.index))

parsers["AmoebaWcaDispersionForce"] = AmoebaWcaDispersionGenerator.parseElement

#=============================================================================================

class AmoebaGeneralizedKirkwoodGenerator:

    """A AmoebaGeneralizedKirkwoodGenerator constructs a AmoebaGeneralizedKirkwoodForce."""
    
    #=========================================================================================

Peter Eastman's avatar
Peter Eastman committed
3564
3565
3566
3567
3568
3569
3570
3571
3572
3573
3574
3575
3576
3577
3578
3579
3580
3581
3582
3583
3584
3585
3586
3587
3588
3589
3590
3591
3592
3593
3594
3595
3596
3597
3598
3599
3600
3601
    def __init__(self, forceField, solventDielectric, soluteDielectric, includeCavityTerm, probeRadius, surfaceAreaFactor):

        self.forceField = forceField
        self.solventDielectric = solventDielectric
        self.soluteDielectric = soluteDielectric
        self.includeCavityTerm = includeCavityTerm
        self.probeRadius = probeRadius
        self.surfaceAreaFactor = surfaceAreaFactor
        self.hasBeenCalled = 0

        self.radiusTypeMap = {}
        self.radiusTypeMap['Bondi'] = {}
        bondiMap = self.radiusTypeMap['Bondi'] 
        rscale = 1.03

        bondiMap[0] = 0.00
        bondiMap[1] = 0.12*rscale
        bondiMap[2] = 0.14*rscale
        bondiMap[5] = 0.18*rscale

        bondiMap[6] = 0.170*rscale
        bondiMap[7] = 0.155*rscale
        bondiMap[8] = 0.152*rscale
        bondiMap[9] = 0.147*rscale

        bondiMap[10] = 0.154*rscale
        bondiMap[14] = 0.210*rscale
        bondiMap[15] = 0.180*rscale
        bondiMap[16] = 0.180*rscale

        bondiMap[17] = 0.175 *rscale
        bondiMap[18] = 0.188*rscale
        bondiMap[34] = 0.190*rscale
        bondiMap[35] = 0.185*rscale

        bondiMap[36] = 0.202*rscale
        bondiMap[53] = 0.198*rscale
        bondiMap[54] = 0.216*rscale
3602
3603
3604

    #=========================================================================================

Peter Eastman's avatar
Peter Eastman committed
3605
    def getObcShct(self, data, atomIndex):
3606

Peter Eastman's avatar
Peter Eastman committed
3607
        atom = data.atoms[atomIndex]
3608
        atomicNumber = atom.element.atomic_number
Peter Eastman's avatar
Peter Eastman committed
3609
        shct = -1.0
3610
3611
3612

        # shct
 
Peter Eastman's avatar
Peter Eastman committed
3613
        if (atomicNumber == 1):                 # H(1)
3614
            shct = 0.85       
Peter Eastman's avatar
Peter Eastman committed
3615
        elif (atomicNumber == 6):               # C(6)
3616
            shct = 0.72         
Peter Eastman's avatar
Peter Eastman committed
3617
        elif (atomicNumber == 7):               # N(7)
3618
            shct = 0.79        
Peter Eastman's avatar
Peter Eastman committed
3619
        elif (atomicNumber == 8):               # O(8)
3620
            shct = 0.85       
Peter Eastman's avatar
Peter Eastman committed
3621
        elif (atomicNumber == 9):               # F(9)
3622
            shct = 0.88    
Peter Eastman's avatar
Peter Eastman committed
3623
        elif (atomicNumber == 15):              # P(15)              
3624
            shct = 0.86 
Peter Eastman's avatar
Peter Eastman committed
3625
        elif (atomicNumber == 16):              # S(16)
3626
            shct = 0.96
Peter Eastman's avatar
Peter Eastman committed
3627
        elif (atomicNumber == 26):              # Fe(26)
3628
3629
            shct = 0.88

Peter Eastman's avatar
Peter Eastman committed
3630
        if (shct < 0.0): 
3631
            shct = 0.80
Peter Eastman's avatar
Peter Eastman committed
3632
            print "getObcShct: Warning no GK overlap scale factor for atom %s of %s %d using default value=%f" % (atom.name, atom.residue.name, atom.residue.index, shct)
3633
3634
3635
3636
3637
 
        return shct 

    #=========================================================================================

Peter Eastman's avatar
Peter Eastman committed
3638
    def getAmoebaTypeRadius(self, data, bondedAtomIndices, atomIndex):
3639

Peter Eastman's avatar
Peter Eastman committed
3640
        atom = data.atoms[atomIndex]
3641
        atomicNumber = atom.element.atomic_number
Peter Eastman's avatar
Peter Eastman committed
3642
        radius = -1.0
3643

Peter Eastman's avatar
Peter Eastman committed
3644
        if (atomicNumber == 1):                  # H(1)
3645
 
Peter Eastman's avatar
Peter Eastman committed
3646
            radius = 0.132
3647
 
Peter Eastman's avatar
Peter Eastman committed
3648
3649
3650
            if (len(bondedAtomIndices) < 1):
                 outputString = "AmoebaGeneralizedKirkwoodGenerator: error getting atom bonded to %s of %s %d " % (atom.name, atom.residue.name, atom.residue.index)
                 raise ValueError(outputString) 
3651
3652
 
            for bondedAtomIndex in bondedAtomIndices:
Peter Eastman's avatar
Peter Eastman committed
3653
                bondedAtomAtomicNumber = data.atoms[bondedAtomIndex].element.atomic_number
3654

Peter Eastman's avatar
Peter Eastman committed
3655
            if (bondedAtomAtomicNumber == 7):
3656
                radius = 0.11
Peter Eastman's avatar
Peter Eastman committed
3657
            if (bondedAtomAtomicNumber == 8):
3658
3659
                radius = 0.105
 
Peter Eastman's avatar
Peter Eastman committed
3660
        elif (atomicNumber == 3):               # Li(3)
3661
            radius = 0.15
Peter Eastman's avatar
Peter Eastman committed
3662
        elif (atomicNumber == 6):               # C(6)
3663
3664
            
            radius = 0.20
Peter Eastman's avatar
Peter Eastman committed
3665
            if (len(bondedAtomIndices) == 3):
3666
3667
                radius = 0.205

Peter Eastman's avatar
Peter Eastman committed
3668
            elif (len(bondedAtomIndices) == 4):
3669
3670
                for bondedAtomIndex in bondedAtomIndices:
                   bondedAtomAtomicNumber = data.atoms[bondedAtomIndex].element.atomic_number
Peter Eastman's avatar
Peter Eastman committed
3671
                   if (bondedAtomAtomicNumber == 7 or bondedAtomAtomicNumber == 8):
3672
3673
                       radius = 0.175

Peter Eastman's avatar
Peter Eastman committed
3674
        elif (atomicNumber == 7):               # N(7)
3675
            radius = 0.16
Peter Eastman's avatar
Peter Eastman committed
3676
        elif (atomicNumber == 8):               # O(8)
3677
            radius = 0.155
Peter Eastman's avatar
Peter Eastman committed
3678
            if (len(bondedAtomIndices) == 2):
3679
                radius = 0.145
Peter Eastman's avatar
Peter Eastman committed
3680
        elif (atomicNumber == 9):               # F(9)
3681
            radius = 0.154
Peter Eastman's avatar
Peter Eastman committed
3682
        elif (atomicNumber == 10):              
3683
            radius = 0.146
Peter Eastman's avatar
Peter Eastman committed
3684
        elif (atomicNumber == 11):              
3685
            radius = 0.209
Peter Eastman's avatar
Peter Eastman committed
3686
        elif (atomicNumber == 12):              
3687
            radius = 0.179
Peter Eastman's avatar
Peter Eastman committed
3688
        elif (atomicNumber == 14):              
3689
            radius = 0.189
Peter Eastman's avatar
Peter Eastman committed
3690
        elif (atomicNumber == 15):              # P(15)              
3691
            radius = 0.196
Peter Eastman's avatar
Peter Eastman committed
3692
        elif (atomicNumber == 16):              # S(16)
3693
            radius = 0.186
Peter Eastman's avatar
Peter Eastman committed
3694
        elif (atomicNumber == 17):              
3695
            radius = 0.182
Peter Eastman's avatar
Peter Eastman committed
3696
        elif (atomicNumber == 18):              
3697
            radius = 0.179
Peter Eastman's avatar
Peter Eastman committed
3698
        elif (atomicNumber == 19):              
3699
            radius = 0.223
Peter Eastman's avatar
Peter Eastman committed
3700
        elif (atomicNumber == 20):              
3701
            radius = 0.191
Peter Eastman's avatar
Peter Eastman committed
3702
        elif (atomicNumber == 35):         
3703
            radius = 2.00
Peter Eastman's avatar
Peter Eastman committed
3704
        elif (atomicNumber == 36):   
3705
            radius = 0.190
Peter Eastman's avatar
Peter Eastman committed
3706
        elif (atomicNumber == 37):              
3707
            radius = 0.226
Peter Eastman's avatar
Peter Eastman committed
3708
        elif (atomicNumber == 53):              
3709
            radius = 0.237
Peter Eastman's avatar
Peter Eastman committed
3710
        elif (atomicNumber == 54):              
3711
            radius = 0.207
Peter Eastman's avatar
Peter Eastman committed
3712
        elif (atomicNumber == 55):              
3713
            radius = 0.263
Peter Eastman's avatar
Peter Eastman committed
3714
        elif (atomicNumber == 56):         
3715
3716
            radius = 0.230

Peter Eastman's avatar
Peter Eastman committed
3717
        if (radius < 0.0): 
3718
            radius = 2.0
Peter Eastman's avatar
Peter Eastman committed
3719
            print "Warning no GK radius for atom %s of %s %d using default value=%f" % (atom.name, atom.residue.name, atom.residue.index, radius)
3720
3721
3722
3723
3724
 
        return radius

    #=========================================================================================

Peter Eastman's avatar
Peter Eastman committed
3725
    def getBondiTypeRadius(self, data, bondedAtomIndices, atomIndex):
3726

Peter Eastman's avatar
Peter Eastman committed
3727
3728
        bondiMap = self.radiusTypeMap['Bondi'] 
        atom = data.atoms[atomIndex]
3729
        atomicNumber = atom.element.atomic_number
Peter Eastman's avatar
Peter Eastman committed
3730
        if (atomicNumber in bondiMap): 
3731
3732
3733
            radius = bondiMap[atomicNumber]
        else:
            radius = 0.206
Peter Eastman's avatar
Peter Eastman committed
3734
            print "Warning no Bondi radius for atom %s of %s %d using default value=%f" % (atom.name, atom.residue.name, atom.residue.index, radius)
3735
3736
3737
3738
3739
3740
3741
3742
3743
3744
3745
3746
 
        return radius

    #=========================================================================================

    @staticmethod
    def parseElement(element, forceField):

        #  <AmoebaGeneralizedKirkwoodForce solventDielectric="78.3" soluteDielectric="1.0" includeCavityTerm="1" probeRadius="0.14" surfaceAreaFactor="-170.351730663">
        #   <GeneralizedKirkwood type="1" charge="-0.22620" shct="0.79"  />
        #   <GeneralizedKirkwood type="2" charge="-0.15245" shct="0.72"  />
        
Peter Eastman's avatar
Peter Eastman committed
3747
        generator = AmoebaGeneralizedKirkwoodGenerator(forceField, element.attrib['solventDielectric'], element.attrib['soluteDielectric'],
3748
                                                        element.attrib['includeCavityTerm'], 
Peter Eastman's avatar
Peter Eastman committed
3749
                                                        element.attrib['probeRadius'], element.attrib['surfaceAreaFactor']) 
3750
3751
3752
3753
        forceField._forces.append(generator)

    #=========================================================================================
    
Peter Eastman's avatar
Peter Eastman committed
3754
    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
3755

Peter Eastman's avatar
Peter Eastman committed
3756
        if (self.hasBeenCalled ):
3757
3758
3759
3760
3761
3762
3763
            return

        verbose = 0
      
        # check if AmoebaMultipoleForce exists since charges needed
        # if it has not been created, raise an error

Peter Eastman's avatar
Peter Eastman committed
3764
        existing = [sys.getForce(i) for i in range(sys.getNumForces())]
3765
        amoebaMultipoleForceList = [f for f in existing if type(f) == mm.AmoebaMultipoleForce]
Peter Eastman's avatar
Peter Eastman committed
3766
        if (len(amoebaMultipoleForceList) > 0):
3767
3768
3769
3770
3771
            amoebaMultipoleForce = amoebaMultipoleForceList[0]
        else:
            # call AmoebaMultipoleForceGenerator.createForce() to ensure charges have been set

            for force in self.forceField._forces:
Peter Eastman's avatar
Peter Eastman committed
3772
3773
                if (force.__class__.__name__ == 'AmoebaMultipoleGenerator'): 
                    force.createForce(sys, data, nonbondedMethod, nonbondedCutoff, args)
3774
3775
3776
3777
3778
3779
3780
3781
3782

        # get or create force depending on whether it has already been added to the system

        existing = [f for f in existing if type(f) == mm.AmoebaGeneralizedKirkwoodForce]
        if len(existing) == 0:

            force = mm.AmoebaGeneralizedKirkwoodForce()
            sys.addForce(force)
 
Peter Eastman's avatar
Peter Eastman committed
3783
3784
            if ('solventDielectric' in args):
                force.setSolventDielectric(float(args['solventDielectric']))
3785
            else:
Peter Eastman's avatar
Peter Eastman committed
3786
                force.setSolventDielectric(   float(self.solventDielectric))
3787

Peter Eastman's avatar
Peter Eastman committed
3788
3789
            if ('soluteDielectric' in args):
                force.setSoluteDielectric(float(args['soluteDielectric']))
3790
            else:
Peter Eastman's avatar
Peter Eastman committed
3791
                force.setSoluteDielectric(    float(self.soluteDielectric))
3792

Peter Eastman's avatar
Peter Eastman committed
3793
3794
            if ('includeCavityTerm' in args):
                force.setIncludeCavityTerm(int(args['includeCavityTerm']))
3795
            else:
Peter Eastman's avatar
Peter Eastman committed
3796
               force.setIncludeCavityTerm(   int(self.includeCavityTerm))
3797
3798
3799
3800
3801
3802
3803
3804
3805
3806

        else:
            print "AmoebaGeneralizedKirkwoodForce exists"
            force = existing[0]

        self.hasBeenCalled = 1

        # add particles to force
        # throw error if particle type not available 

Peter Eastman's avatar
Peter Eastman committed
3807
3808
        force.setProbeRadius(         float(self.probeRadius))
        force.setSurfaceAreaFactor(   float(self.surfaceAreaFactor))
3809
3810
3811
3812
3813

        # 1-2

        bonded12ParticleSets = []
        for i in range(len(data.atoms)):
Peter Eastman's avatar
Peter Eastman committed
3814
3815
3816
            bonded12ParticleSet = AmoebaVdwGenerator.getBondedParticleSet(i, data)
            bonded12ParticleSet = set(sorted(bonded12ParticleSet))
            bonded12ParticleSets.append(bonded12ParticleSet)
3817
3818

        radiusType = 'Bondi'
Peter Eastman's avatar
Peter Eastman committed
3819
3820
3821
3822
        for atomIndex in range(0, amoebaMultipoleForce.getNumMultipoles()):
            multipoleParameters = amoebaMultipoleForce.getMultipoleParameters(atomIndex)
            if (radiusType == 'Amoeba'):
                radius = self.getAmoebaTypeRadius(data, bonded12ParticleSets[atomIndex], atomIndex)
3823
            else:
Peter Eastman's avatar
Peter Eastman committed
3824
3825
3826
3827
3828
                radius = self.getBondiTypeRadius(data, bonded12ParticleSets[atomIndex], atomIndex)
            shct = self.getObcShct(data, atomIndex)
            force.addParticle(multipoleParameters[0], radius, shct)
            if (verbose):
               print "GeneralizedKirkwood %5d %15.7e %15.7e" % (atomIndex, multipoleParameters[0], radius, shct)
3829
3830
3831
3832
3833
3834
3835
3836
3837
3838
3839
3840
3841
3842
3843

parsers["AmoebaGeneralizedKirkwoodForce"] = AmoebaGeneralizedKirkwoodGenerator.parseElement

#=============================================================================================

class AmoebaUreyBradleyGenerator:

    #=============================================================================================

    """An AmoebaUreyBradleyGenerator constructs a AmoebaUreyBradleyForce."""

    #=============================================================================================
    
    def __init__(self, cubic, quartic):

Peter Eastman's avatar
Peter Eastman committed
3844
3845
        self.cubic = cubic
        self.quartic = quartic
3846

Peter Eastman's avatar
Peter Eastman committed
3847
3848
3849
        self.types1 = []
        self.types2 = []
        self.types3 = []
3850

Peter Eastman's avatar
Peter Eastman committed
3851
3852
        self.length = []
        self.k = []
3853
3854
3855
3856
3857
3858
3859
3860
3861
3862
3863

        self.hasBeenCalled = 0
    
    #=============================================================================================

    @staticmethod
    def parseElement(element, forceField):

        #  <AmoebaUreyBradleyForce cubic="0.0" quartic="0.0"  >
        #   <UreyBradley class1="74" class2="73" class3="74" k="16003.8" d="0.15537" /> 

Peter Eastman's avatar
Peter Eastman committed
3864
        generator = AmoebaUreyBradleyGenerator(float(element.attrib['cubic']), float(element.attrib['quartic']))
3865
3866
3867
3868
3869
3870
3871
3872
3873
3874
3875
3876
3877
3878
        forceField._forces.append(generator)
        for bond in element.findall('UreyBradley'):
            types = forceField._findAtomTypes(bond, 3)
            if types is not None:

                generator.types1.append(types[0])
                generator.types2.append(types[1])
                generator.types3.append(types[2])

                generator.length.append(float(bond.attrib['d']))
                generator.k.append(float(bond.attrib['k']))

            else:
                outputString = "AmoebaUreyBradleyGenerator: error getting types: %s %s %s" % (
Peter Eastman's avatar
Peter Eastman committed
3879
3880
                                    bond.attrib['class1'], bond.attrib['class2'], bond.attrib['class3'])
                raise ValueError(outputString) 
3881
3882
3883
    
    #=============================================================================================

Peter Eastman's avatar
Peter Eastman committed
3884
    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
3885

Peter Eastman's avatar
Peter Eastman committed
3886
        verbose = 0
3887

Peter Eastman's avatar
Peter Eastman committed
3888
        if (self.hasBeenCalled):
3889
3890
3891
3892
             return
 
        self.hasBeenCalled = 1

Peter Eastman's avatar
Peter Eastman committed
3893
3894
        existing = [sys.getForce(i) for i in range(sys.getNumForces())]
        existing = [f for f in existing if type(f) == mm.AmoebaUreyBradleyForce]
3895
3896

        if len(existing) == 0:
Peter Eastman's avatar
Peter Eastman committed
3897
            force = mm.AmoebaUreyBradleyForce()
3898
3899
3900
3901
            sys.addForce(force)
        else:
            force = existing[0]

Peter Eastman's avatar
Peter Eastman committed
3902
3903
        force.setAmoebaGlobalUreyBradleyCubic(self.cubic)
        force.setAmoebaGlobalUreyBradleyQuartic(self.quartic)
3904
3905

        for (angle, isConstrained) in zip(data.angles, data.isAngleConstrained):
Peter Eastman's avatar
Peter Eastman committed
3906
            if (isConstrained):
3907
3908
3909
3910
3911
3912
3913
3914
                continue
            type1 = data.atomType[data.atoms[angle[0]]]
            type2 = data.atomType[data.atoms[angle[1]]]
            type3 = data.atomType[data.atoms[angle[2]]]
            for i in range(len(self.types1)):
                types1 = self.types1[i]
                types2 = self.types2[i]
                types3 = self.types3[i]
Peter Eastman's avatar
Peter Eastman committed
3915
                if ((type1 in types1 and type2 in types2 and type3 in types3) or (type3 in types1 and type2 in types2 and type1 in types3)):
3916

Peter Eastman's avatar
Peter Eastman committed
3917
                    if (verbose):
3918
3919
3920
3921
3922
3923
3924
3925
                        print "AmoebaUreyBradleyGenerator %5d %5d %5d [%5s %5s %5s] %15.6f %15.6f" % (angle[0], angle[1], angle[2], type1, type2, type3, self.length[i], self.k[i])

                    force.addUreyBradley(angle[0], angle[2], self.length[i], self.k[i])
                    break

parsers["AmoebaUreyBradleyForce"] = AmoebaUreyBradleyGenerator.parseElement

#=============================================================================================