@@ -407,142 +407,3 @@ def create_system(self, topology, molecules=None):
407407 system = self .postprocess_system (system )
408408
409409 return system
410-
411-
412- ################################################################################
413- # Dummy system generator
414- ################################################################################
415-
416-
417- class DummySystemGenerator (SystemGenerator ):
418- """
419- Dummy force field that can add basic parameters to any system for testing purposes.
420-
421- * All particles are assigned carbon mass
422- * All particles interact with a repulsive potential
423- * All bonds have equilibrium length 1 A
424- * All angles have equilibrium angle dependent on number of substituents of central atom
425- * 2, 3 bonds: 120 degrees
426- * 4 bonds: 109.8 degrees
427- * 5 or more bonds: 90 degrees
428- * Torsions are added with periodicity 3, but no barrier height
429-
430- """
431-
432- def create_system (self , topology , ** kwargs ):
433- """
434- Create a System object with simple parameters from the provided Topology
435-
436- Any kwargs are ignored.
437-
438- Parameters
439- ----------
440- topology : openff.toolkit.topology.Topology
441- The Topology to be parameterized
442-
443- Returns
444- -------
445- system : openmm.System
446- The System object
447-
448- """
449- # TODO: Allow periodicity to be determined from topology
450-
451- import openmm
452- from openmm import unit
453- from openmmtools .constants import kB
454-
455- kT = kB * 300 * unit .kelvin # hard-coded temperature for setting energy scales
456-
457- # Create a System
458- system = openmm .System ()
459-
460- # Add particles
461- mass = 12.0 * unit .amu
462- for atom in topology .atoms :
463- system .addParticle (mass )
464-
465- # Add simple repulsive interactions
466- # TODO: Use softcore repulsive interaction; Gaussian times switch?
467- nonbonded = openmm .CustomNonbondedForce ("100/(r/0.1)^4" )
468- nonbonded .setNonbondedMethod (openmm .CustomNonbondedForce .CutoffNonPeriodic )
469- nonbonded .setCutoffDistance (1 * unit .nanometer )
470- system .addForce (nonbonded )
471- for atom in topology .atoms :
472- nonbonded .addParticle ([])
473-
474- # Build a list of which atom indices are bonded to each atom
475- bondedToAtom = []
476- for atom in topology .atoms ():
477- bondedToAtom .append (set ())
478- for atom1 , atom2 in topology .bonds ():
479- bondedToAtom [atom1 .index ].add (atom2 .index )
480- bondedToAtom [atom2 .index ].add (atom1 .index )
481- return bondedToAtom
482-
483- # Add bonds
484- bond_force = openmm .HarmonicBondForce ()
485- r0 = 1.0 * unit .angstroms
486- sigma_r = 0.1 * unit .angstroms
487- Kr = kT / sigma_r ** 2
488- for atom1 , atom2 in topology .bonds ():
489- bond_force .addBond (atom1 .index , atom2 .index , r0 , Kr )
490- system .addForce (bond_force )
491-
492- # Add angles
493- uniqueAngles = set ()
494- for bond in topology .bonds ():
495- for atom in bondedToAtom [bond .atom1 ]:
496- if atom != bond .atom2 :
497- if atom < bond .atom2 :
498- uniqueAngles .add ((atom , bond .atom1 , bond .atom2 ))
499- else :
500- uniqueAngles .add ((bond .atom2 , bond .atom1 , atom ))
501- for atom in bondedToAtom [bond .atom2 ]:
502- if atom != bond .atom1 :
503- if atom > bond .atom1 :
504- uniqueAngles .add ((bond .atom1 , bond .atom2 , atom ))
505- else :
506- uniqueAngles .add ((atom , bond .atom2 , bond .atom1 ))
507- angles = sorted (list (uniqueAngles ))
508- theta0 = 109.5 * unit .degrees # TODO: Adapt based on number of bonds to each atom?
509- sigma_theta = 10 * unit .degrees
510- Ktheta = kT / sigma_theta ** 2
511- angle_force = openmm .HarmonicAngleForce ()
512- for atom1 , atom2 , atom3 in angles :
513- angles .addAngle (atom1 .index , atom2 .index , atom3 .index , theta0 , Ktheta )
514- system .addForce (angle_force )
515-
516- # Make a list of all unique proper torsions
517- uniquePropers = set ()
518- for angle in angles :
519- for atom in bondedToAtom [angle [0 ]]:
520- if atom not in angle :
521- if atom < angle [2 ]:
522- uniquePropers .add ((atom , angle [0 ], angle [1 ], angle [2 ]))
523- else :
524- uniquePropers .add ((angle [2 ], angle [1 ], angle [0 ], atom ))
525- for atom in bondedToAtom [angle [2 ]]:
526- if atom not in angle :
527- if atom > angle [0 ]:
528- uniquePropers .add ((angle [0 ], angle [1 ], angle [2 ], atom ))
529- else :
530- uniquePropers .add ((atom , angle [2 ], angle [1 ], angle [0 ]))
531- propers = sorted (list (uniquePropers ))
532- torsion_force = openmm .PeriodicTorsionForce ()
533- periodicity = 3
534- phase = 0.0 * unit .degrees
535- Kphi = 0.0 * kT
536- for atom1 , atom2 , atom3 , atom4 in propers :
537- torsion_force .add_torsion (
538- atom1 .index ,
539- atom2 .index ,
540- atom3 .index ,
541- atom4 .index ,
542- periodicity ,
543- phase ,
544- Kphi ,
545- )
546- system .addForce (torsion_force )
547-
548- return system
0 commit comments