"""
This module contains a chamber prmtop class that will read in all
parameters and allow users to manipulate that data and write a new
prmtop object.
Copyright (C) 2010 - 2015 Jason Swails
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330,
Boston, MA 02111-1307, USA.
"""
from __future__ import absolute_import, division, print_function
import copy as _copy
import warnings
from math import pi, sqrt
from ..constants import DEG_TO_RAD, IFBOX, NATOM, NATYP, NTYPES, RAD_TO_DEG, SMALL, TINY
from ..exceptions import AmberError, AmberWarning
from ..topologyobjects import BondType, ExtraPoint, Improper, ImproperType, UreyBradley
from ..utils.six.moves import range, zip
from ._amberparm import AmberParm
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[docs]class ChamberParm(AmberParm):
"""Chamber Topology (parm7 format) class.
Gives low, and some high, level access to topology data or
interact with some of the high-level classes comprising the system
topology and parameters. The ChamberParm class uses the same
attributes that the AmberParm class uses, and only the ones unique
to ChamberParm will be shown below.
Parameters
----------
prm_name : str, optional
If provided, this file is parsed and the data structures will be loaded
from the data in this file
xyz : str or array, optional
If provided, the coordinates and unit cell dimensions from the provided
Amber inpcrd/restart file will be loaded into the molecule, or the
coordinates will be loaded from the coordinate array
box : array, optional
If provided, the unit cell information will be set from the provided
unit cell dimensions (a, b, c, alpha, beta, and gamma, respectively)
Attributes
----------
LJ_14_radius : list(float)
The same as LJ_radius, except specific for 1-4 nonbonded parameters,
which may differ in the CHARMM force field
LJ_14_depth : list(float)
The same as LJ_depth, except specific for 1-4 nonbonded parameters,
which may differ in the CHARMM force field
urey_bradleys : TrackedList(UreyBradley)
List of Urey-Bradley terms between two atoms in a valence angle
impropers : TrackedList(Improper)
List of CHARMM-style improper torsions
cmaps : TrackedList(Cmap)
List of coupled-torsion correction map parameters
urey_bradley_types : TrackedList(UreyBradleyType)
List of parameters defining the Urey-Bradley terms
improper_types : TrackedList(Improper)
List of parameters defining the Improper terms
cmap_types : TrackedList(CmapType)
List of parameters defining the CMAP terms
chamber : bool=True
On ChamberParm instances, this is always True to indicate that it is a
CHAMBER-style topology file
amoeba : bool=False
On ChamberParm instances, this is always False to indicate that it is
not an AMOEBA-style topology file
has_cmap : bool
True if CMAP parameters are present in this system; False otherwise
See Also
--------
:class:`AmberParm`
"""
_cmap_prefix = "CHARMM_"
#===================================================
[docs] def initialize_topology(self, xyz=None, box=None):
"""
Initializes topology data structures, like the list of atoms, bonds,
etc., after the topology file has been read. The following methods are
called:
"""
self.LJ_14_radius = []
self.LJ_14_depth = []
AmberParm.initialize_topology(self, xyz, box)
#===================================================
def _copy_lj_data(self, other):
""" Copies Lennard-Jones lists and dicts from myself to a copy """
super(ChamberParm, self)._copy_lj_data(other)
other.LJ_14_radius = _copy.copy(self.LJ_14_radius)
other.LJ_14_depth = _copy.copy(self.LJ_14_depth)
for atom in other.atoms:
other.LJ_14_radius[atom.nb_idx-1] = atom.atom_type.rmin_14
other.LJ_14_depth[atom.nb_idx-1] = atom.atom_type.epsilon_14
#===================================================
[docs] def load_pointers(self):
"""
Loads the data in POINTERS section into a pointers dictionary with each
key being the pointer name according to http://ambermd.org/formats.html
"""
AmberParm.load_pointers(self)
# Other pointers
nub, nubtypes = self.parm_data['CHARMM_UREY_BRADLEY_COUNT'][:2]
self.pointers['NUB'] = nub
self.pointers['NUBTYPES'] = nubtypes
self.pointers['NIMPHI'] = self.parm_data['CHARMM_NUM_IMPROPERS'][0]
self.pointers['NIMPRTYPES'] = self.parm_data['CHARMM_NUM_IMPR_TYPES'][0]
#===================================================
[docs] def load_structure(self):
"""
Loads all of the topology instance variables. This is necessary if we
actually want to modify the topological layout of our system
(like deleting atoms)
"""
super(ChamberParm, self).load_structure()
self._load_urey_brad_info()
self._load_improper_info()
# All of our angles are urey-bradley types
for angle in self.angles: angle.funct = 5
super(ChamberParm, self).unchange()
#===================================================
[docs] @classmethod
def from_structure(cls, struct, copy=False):
"""
Take a Structure instance and initialize a ChamberParm instance from
that data.
Parameters
----------
struct : Structure
The input structure from which to construct a ChamberParm instance
copy : bool
If True, the input struct is deep-copied to make sure it does not
share any objects with the original ``struct``. Default is False
Returns
-------
inst : :class:`ChamberParm`
The ChamberParm instance derived from the input structure
Notes
-----
Due to the nature of the prmtop file, struct almost *always* returns a
deep copy. The one exception is when struct is already of type
:class:`ChamberParm`, in which case the original object is returned
unless ``copy`` is ``True``.
"""
if isinstance(struct, cls):
if copy:
return _copy.copy(struct)
return struct
if (struct.rb_torsions or struct.trigonal_angles or struct.pi_torsions
or struct.out_of_plane_bends or struct.stretch_bends
or struct.torsion_torsions or struct.multipole_frames):
raise TypeError('ChamberParm does not support all potential terms '
'defined in the input Structure')
inst = struct.copy(cls, split_dihedrals=True)
inst.update_dihedral_exclusions()
inst._add_missing_13_14(ignore_inconsistent_vdw=True)
inst.pointers = {}
inst.LJ_types = {}
nbfixes = inst.atoms.assign_nbidx_from_types()
# Give virtual sites a name that Amber understands
for atom in inst.atoms:
if isinstance(atom, ExtraPoint): atom.type = 'EP'
# Fill the Lennard-Jones arrays/dicts
ntyp = 0
for atom in inst.atoms:
inst.LJ_types[atom.type] = atom.nb_idx
ntyp = max(ntyp, atom.nb_idx)
inst.LJ_radius = [0 for i in range(ntyp)]
inst.LJ_depth = [0 for i in range(ntyp)]
inst.LJ_14_radius = [0 for i in range(ntyp)]
inst.LJ_14_depth = [0 for i in range(ntyp)]
for atom in inst.atoms:
inst.LJ_radius[atom.nb_idx-1] = atom.rmin
inst.LJ_depth[atom.nb_idx-1] = atom.epsilon
inst.LJ_14_radius[atom.nb_idx-1] = atom.rmin_14
inst.LJ_14_depth[atom.nb_idx-1] = atom.epsilon_14
inst._add_standard_flags()
inst.pointers['NATOM'] = len(inst.atoms)
inst.parm_data['POINTERS'][NATOM] = len(inst.atoms)
inst.box = _copy.copy(struct.box)
if struct.box is None:
inst.parm_data['POINTERS'][IFBOX] = 0
inst.pointers['IFBOX'] = 0
elif (abs(struct.box[3] - 90) > TINY or abs(struct.box[4] - 90) > TINY
or abs(struct.box[5] - 90) > TINY):
inst.parm_data['POINTERS'][IFBOX] = 2
inst.pointers['IFBOX'] = 2
inst.parm_data['BOX_DIMENSIONS'] = [struct.box[3]] + list(struct.box[:3])
else:
inst.parm_data['POINTERS'][IFBOX] = 1
inst.pointers['IFBOX'] = 1
inst.parm_data['BOX_DIMENSIONS'] = [90] + list(struct.box[:3])
# pmemd likes to skip torsions with periodicities of 0, which may be
# present as a way to hack entries into the 1-4 pairlist. See
# https://github.com/ParmEd/ParmEd/pull/145 for discussion. The solution
# here is to simply set that periodicity to 1.
for dt in inst.dihedral_types:
if dt.phi_k == 0 and dt.per == 0:
dt.per = 1.0
elif dt.per == 0:
warnings.warn('Periodicity of 0 detected with non-zero force constant. Changing '
'periodicity to 1 and force constant to 0 to ensure 1-4 nonbonded '
'pairs are properly identified. This might cause a shift in the '
'energy, but will leave forces unaffected', AmberWarning)
dt.phi_k = 0.0
dt.per = 1.0
inst.remake_parm()
inst._set_nonbonded_tables(nbfixes)
del inst.adjusts[:]
inst.parm_data['FORCE_FIELD_TYPE'] = fftype = []
fftype.extend([1, 'CHARMM force field: No FF information parsed...'])
return inst
#===================================================
[docs] def remake_parm(self):
"""
Fills :attr:`parm_data` from the data in the parameter and topology
arrays (e.g., :attr:`atoms`, :attr:`bonds`, :attr:`bond_types`, ...)
"""
# Get rid of terms containing deleted atoms and empty residues
self.prune_empty_terms()
self.residues.prune()
self.rediscover_molecules()
# Transfer information from the topology lists
self._xfer_atom_info()
self._xfer_residue_info()
self._xfer_bond_info()
self._xfer_angle_info()
self._xfer_dihedral_info()
self._xfer_urey_bradley_properties()
self._xfer_improper_properties()
self._xfer_cmap_properties()
# Load the pointers dict
self.load_pointers()
# Mark atom list as unchanged
super(ChamberParm, self).unchange()
#===================================================
[docs] def fill_LJ(self):
"""
Fills the LJ_radius, LJ_depth arrays and LJ_types dictionary with data
from LENNARD_JONES_ACOEF and LENNARD_JONES_BCOEF sections of the prmtop
files, by undoing the canonical combining rules.
"""
AmberParm.fill_LJ(self)
acoef = self.parm_data['LENNARD_JONES_14_ACOEF']
bcoef = self.parm_data['LENNARD_JONES_14_BCOEF']
ntypes = self.pointers['NTYPES']
self.LJ_14_radius = [] # empty LJ_radii so it can be re-filled
self.LJ_14_depth = [] # empty LJ_depths so it can be re-filled
one_sixth = 1.0 / 6.0 # we need to raise some numbers to the 1/6th power
for i in range(ntypes):
lj_index = self.parm_data["NONBONDED_PARM_INDEX"][ntypes*i+i] - 1
if acoef[lj_index] < 1.0e-6:
self.LJ_14_radius.append(0)
self.LJ_14_depth.append(0)
else:
factor = 2 * acoef[lj_index] / bcoef[lj_index]
self.LJ_14_radius.append(pow(factor, one_sixth) * 0.5)
self.LJ_14_depth.append(bcoef[lj_index] / 2 / factor)
#===================================================
[docs] def recalculate_LJ(self):
"""
Takes the values of the LJ_radius and LJ_depth arrays and recalculates
the LENNARD_JONES_A/BCOEF topology sections from the canonical
combining rules.
"""
AmberParm.recalculate_LJ(self)
ntypes = self.pointers['NTYPES']
acoef = self.parm_data['LENNARD_JONES_14_ACOEF']
bcoef = self.parm_data['LENNARD_JONES_14_BCOEF']
for i in range(ntypes):
for j in range(i, ntypes):
index = self.parm_data['NONBONDED_PARM_INDEX'][ntypes*i+j] - 1
rij = self.LJ_14_radius[i] + self.LJ_14_radius[j]
wdij = sqrt(self.LJ_14_depth[i] * self.LJ_14_depth[j])
acoef[index] = wdij * rij ** 12
bcoef[index] = 2 * wdij * rij**6
#===================================================
@property
def chamber(self):
return True
@property
def amoeba(self):
return False
#=========== PRIVATE INSTANCE METHODS ============
def _load_urey_brad_info(self):
""" Loads the Urey-Bradley types and array """
del self.urey_bradleys[:]
del self.urey_bradley_types[:]
for k, req in zip(self.parm_data['CHARMM_UREY_BRADLEY_FORCE_CONSTANT'],
self.parm_data['CHARMM_UREY_BRADLEY_EQUIL_VALUE']):
self.urey_bradley_types.append(
BondType(k, req, self.urey_bradley_types)
)
it = iter(self.parm_data['CHARMM_UREY_BRADLEY'])
for i, j, k in zip(it, it, it):
self.urey_bradleys.append(
UreyBradley(self.atoms[i-1], self.atoms[j-1],
self.urey_bradley_types[k-1])
)
#===================================================
def _load_improper_info(self):
""" Loads the CHARMM Improper types and array """
del self.impropers[:]
del self.improper_types[:]
for k, eq in zip(self.parm_data['CHARMM_IMPROPER_FORCE_CONSTANT'],
self.parm_data['CHARMM_IMPROPER_PHASE']):
# Previous versions of ParmEd stored improper phases as degrees,
# whereas it should really be stored in radians. So do a simple
# heuristic check to see if a conversion is necessary so we support
# all versions.
eq = eq * RAD_TO_DEG if abs(eq) <= 2*pi else eq
self.improper_types.append(
ImproperType(k, eq, self.improper_types)
)
it = iter(self.parm_data['CHARMM_IMPROPERS'])
for i, j, k, l, m in zip(it, it, it, it, it):
self.impropers.append(
Improper(self.atoms[i-1], self.atoms[j-1], self.atoms[k-1],
self.atoms[l-1], self.improper_types[m-1])
)
# Make sure that if we have a comment in the CHARMM impropers, we fix it
# to say the units are in radians
for i in range(len(self.parm_comments.get('CHARMM_IMPROPER_PHASE', []))):
comment = self.parm_comments['CHARMM_IMPROPER_PHASE'][i]
if 'degrees' in comment:
self.parm_comments['CHARMM_IMPROPER_PHASE'][i] = \
comment.replace('degrees', 'radians')
#===================================================
def _check_section_lengths(self):
""" Make sure each section has the necessary number of entries """
super(ChamberParm, self)._check_section_lengths()
def check_length(key, length, required=True):
if not required and not key in self.parm_data: return
if len(self.parm_data[key]) != length:
raise AmberError('FLAG %s has %d elements; expected %d' %
(key, len(self.parm_data[key]), length))
check_length('CHARMM_UREY_BRADLEY_COUNT', 2)
check_length('CHARMM_UREY_BRADLEY', self.pointers['NUB']*3)
check_length('CHARMM_UREY_BRADLEY_FORCE_CONSTANT',
self.pointers['NUBTYPES'])
check_length('CHARMM_UREY_BRADLEY_EQUIL_VALUE',
self.pointers['NUBTYPES'])
check_length('CHARMM_NUM_IMPROPERS', 1)
check_length('CHARMM_IMPROPERS', self.pointers['NIMPHI']*5)
check_length('CHARMM_NUM_IMPR_TYPES', 1)
check_length('CHARMM_IMPROPER_FORCE_CONSTANT',
self.pointers['NIMPRTYPES'])
check_length('CHARMM_IMPROPER_PHASE', self.pointers['NIMPRTYPES'])
ntypes = self.pointers['NTYPES']
check_length('LENNARD_JONES_14_ACOEF', ntypes*(ntypes+1)//2)
check_length('LENNARD_JONES_14_BCOEF', ntypes*(ntypes+1)//2)
#===================================================
def _xfer_urey_bradley_properties(self):
"""
Sets the various topology file section data from the Urey-Bradley arrays
"""
data = self.parm_data
for urey_type in self.urey_bradley_types:
urey_type.used = False
for urey in self.urey_bradleys:
urey.type.used = True
self.urey_bradley_types.prune_unused()
data['CHARMM_UREY_BRADLEY_FORCE_CONSTANT'] = \
[type.k for type in self.urey_bradley_types]
data['CHARMM_UREY_BRADLEY_EQUIL_VALUE'] = \
[type.req for type in self.urey_bradley_types]
data['CHARMM_UREY_BRADLEY'] = bond_array = []
for urey in self.urey_bradleys:
bond_array.extend([urey.atom1.idx+1, urey.atom2.idx+1,
urey.type.idx+1])
nub = len(self.urey_bradleys)
nubt = len(self.urey_bradley_types)
data['CHARMM_UREY_BRADLEY_COUNT'] = [nub, nubt]
self.pointers['NUB'] = nub
self.pointers['NUBTYPES'] = nubt
#===================================================
def _xfer_improper_properties(self):
""" Sets the topology file section data from the improper arrays """
data = self.parm_data
for improper_type in self.improper_types:
improper_type.used = False
for improper in self.impropers:
improper.type.used = True
self.improper_types.prune_unused()
data['CHARMM_IMPROPER_FORCE_CONSTANT'] = \
[type.psi_k for type in self.improper_types]
data['CHARMM_IMPROPER_PHASE'] = \
[type.psi_eq*DEG_TO_RAD for type in self.improper_types]
data['CHARMM_IMPROPERS'] = improper_array = []
for imp in self.impropers:
improper_array.extend([imp.atom1.idx+1, imp.atom2.idx+1,
imp.atom3.idx+1, imp.atom4.idx+1,
imp.type.idx+1])
data['CHARMM_NUM_IMPROPERS'] = [len(self.impropers)]
data['CHARMM_NUM_IMPR_TYPES'] = [len(self.improper_types)]
self.pointers['NIMPHI'] = len(improper_array)
self.pointers['NIMPRTYPES'] = len(self.improper_types)
#===================================================
def _add_standard_flags(self):
""" Adds all of the standard flags to the parm_data array """
self.set_version()
self.add_flag('CTITLE', '20a4', num_items=0)
self.add_flag('POINTERS', '10I8', num_items=31)
self.add_flag('FORCE_FIELD_TYPE', 'i2,a78', num_items=0)
self.add_flag('ATOM_NAME', '20a4', num_items=0)
self.add_flag('CHARGE', '3E24.16', num_items=0,
comments=['Atomic charge multiplied by sqrt(332.0716D0) (CCELEC)'])
self.add_flag('ATOMIC_NUMBER', '10I8', num_items=0)
self.add_flag('MASS', '5E16.8', num_items=0)
self.add_flag('ATOM_TYPE_INDEX', '10I8', num_items=0)
self.add_flag('NUMBER_EXCLUDED_ATOMS', '10I8', num_items=0)
self.add_flag('NONBONDED_PARM_INDEX', '10I8', num_items=0)
self.add_flag('RESIDUE_LABEL', '20a4', num_items=0)
self.add_flag('RESIDUE_POINTER', '10I8', num_items=0)
self.add_flag('BOND_FORCE_CONSTANT', '5E16.8', num_items=0)
self.add_flag('BOND_EQUIL_VALUE', '5E16.8', num_items=0)
self.add_flag('ANGLE_FORCE_CONSTANT', '5E16.8', num_items=0)
self.add_flag('ANGLE_EQUIL_VALUE', '3E25.17', num_items=0)
self.add_flag('CHARMM_UREY_BRADLEY_COUNT', '2I8', num_items=2,
comments=['V(ub) = K_ub(r_ik - R_ub)**2',
'Number of Urey Bradley terms and types'])
self.add_flag('CHARMM_UREY_BRADLEY', '10I8', num_items=0,
comments=['List of the two atoms and its parameter index',
'in each UB term: i,k,index'])
self.add_flag('CHARMM_UREY_BRADLEY_FORCE_CONSTANT', '5E16.8', num_items=0,
comments=['K_ub: kcal/mol/A**2'])
self.add_flag('CHARMM_UREY_BRADLEY_EQUIL_VALUE', '5E16.8', num_items=0,
comments=['r_ub: A'])
self.add_flag('DIHEDRAL_FORCE_CONSTANT', '5E16.8', num_items=0)
self.add_flag('DIHEDRAL_PERIODICITY', '5E16.8', num_items=0)
self.add_flag('DIHEDRAL_PHASE', '5E16.8', num_items=0)
self.add_flag('SCEE_SCALE_FACTOR', '5E16.8', num_items=0)
self.add_flag('SCNB_SCALE_FACTOR', '5E16.8', num_items=0)
self.add_flag('CHARMM_NUM_IMPROPERS', '10I8', num_items=0,
comments=['Number of terms contributing to the',
'quadratic four atom improper energy term:',
'V(improper) = K_psi(psi - psi_0)**2'])
self.add_flag('CHARMM_IMPROPERS', '10I8', num_items=0,
comments=['List of the four atoms in each improper term',
'i,j,k,l,index i,j,k,l,index',
'where index is into the following two lists:',
'CHARMM_IMPROPER_{FORCE_CONSTANT,IMPROPER_PHASE}'])
self.add_flag('CHARMM_NUM_IMPR_TYPES', '1I8', num_items=1,
comments=['Number of unique parameters contributing to the',
'quadratic four atom improper energy term'])
self.add_flag('CHARMM_IMPROPER_FORCE_CONSTANT', '5E16.8', num_items=0,
comments=['K_psi: kcal/mole/rad**2'])
self.add_flag('CHARMM_IMPROPER_PHASE', '5E16.8', num_items=0, comments=['psi: radians'])
natyp = self.pointers['NATYP'] = self.parm_data['POINTERS'][NATYP] = 1
self.add_flag('SOLTY', '5E16.8', num_items=natyp)
self.add_flag('LENNARD_JONES_ACOEF', '3E24.16', num_items=0)
self.add_flag('LENNARD_JONES_BCOEF', '3E24.16', num_items=0)
self.add_flag('LENNARD_JONES_14_ACOEF', '3E24.16', num_items=0)
self.add_flag('LENNARD_JONES_14_BCOEF', '3E24.16', num_items=0)
self.add_flag('BONDS_INC_HYDROGEN', '10I8', num_items=0)
self.add_flag('BONDS_WITHOUT_HYDROGEN', '10I8', num_items=0)
self.add_flag('ANGLES_INC_HYDROGEN', '10I8', num_items=0)
self.add_flag('ANGLES_WITHOUT_HYDROGEN', '10I8', num_items=0)
self.add_flag('DIHEDRALS_INC_HYDROGEN', '10I8', num_items=0)
self.add_flag('DIHEDRALS_WITHOUT_HYDROGEN', '10I8', num_items=0)
self.add_flag('EXCLUDED_ATOMS_LIST', '10I8', num_items=0)
self.add_flag('HBOND_ACOEF', '5E16.8', num_items=0)
self.add_flag('HBOND_BCOEF', '5E16.8', num_items=0)
self.add_flag('HBCUT', '5E16.8', num_items=0)
self.add_flag('AMBER_ATOM_TYPE', '20a4', num_items=0)
self.add_flag('TREE_CHAIN_CLASSIFICATION', '20a4', num_items=0)
self.add_flag('JOIN_ARRAY', '10I8', num_items=0)
self.add_flag('IROTAT', '10I8', num_items=0)
if self.has_cmap:
self.add_flag(self._cmap_prefix + 'CMAP_COUNT', '2I8', num_items=2,
comments=['Number of CMAP terms, number of unique CMAP parameters'])
self.add_flag(self._cmap_prefix + 'CMAP_RESOLUTION', '20I4', num_items=0,
comments=['Number of steps along each phi/psi CMAP axis',
'for each CMAP_PARAMETER grid'])
self.add_flag(self._cmap_prefix + 'CMAP_INDEX', '6I8', num_items=0,
comments=['Atom index i,j,k,l,m of the cross term',
'and then pointer to CMAP_PARAMETER_n'])
if self.box is not None:
self.add_flag('SOLVENT_POINTERS', '3I8', num_items=3)
self.add_flag('ATOMS_PER_MOLECULE', '10I8', num_items=0)
self.add_flag('BOX_DIMENSIONS', '5E16.8', num_items=4)
self.add_flag('RADIUS_SET', '1a80', num_items=1)
self.add_flag('RADII', '5E16.8', num_items=0)
self.add_flag('SCREEN', '5E16.8', num_items=0)
self.add_flag('IPOL', '1I8', num_items=1)
#===================================================
def _set_nonbonded_tables(self, nbfixes=None):
""" Sets the tables of Lennard-Jones nonbonded interaction pairs """
from parmed.tools.actions import addLJType
data = self.parm_data
ntypes = data['POINTERS'][NTYPES]
ntypes2 = ntypes * ntypes
# Set up the index lookup tables (not a unique solution)
data['NONBONDED_PARM_INDEX'] = [0 for i in range(ntypes2)]
holder = [0 for i in range(ntypes2)]
idx = 0
for i in range(ntypes):
for j in range(i+1):
idx += 1
holder[ntypes*i+j] = holder[ntypes*j+i] = idx
idx = 0
for i in range(ntypes):
for j in range(ntypes):
data['NONBONDED_PARM_INDEX'][idx] = holder[ntypes*i+j]
idx += 1
nttyp = ntypes * (ntypes + 1) // 2
# Now build the Lennard-Jones arrays
data['LENNARD_JONES_14_ACOEF'] = [0 for i in range(nttyp)]
data['LENNARD_JONES_14_BCOEF'] = [0 for i in range(nttyp)]
data['LENNARD_JONES_ACOEF'] = [0 for i in range(nttyp)]
data['LENNARD_JONES_BCOEF'] = [0 for i in range(nttyp)]
self.recalculate_LJ()
# Now make any NBFIX modifications we had
if nbfixes is not None:
for i, fix in enumerate(nbfixes):
for terms in fix:
j, rmin, eps, rmin14, eps14 = terms
i, j = min(i, j-1), max(i, j-1)
eps = abs(eps)
eps14 = abs(eps14)
idx = data['NONBONDED_PARM_INDEX'][ntypes*i+j] - 1
data['LENNARD_JONES_ACOEF'][idx] = eps * rmin**12
data['LENNARD_JONES_BCOEF'][idx] = 2 * eps * rmin**6
data['LENNARD_JONES_14_ACOEF'][idx] = eps14 * rmin14**12
data['LENNARD_JONES_14_BCOEF'][idx] = 2 * eps14 * rmin14**6
# If we had an explicit set of exceptions, we need to implement all of
# those exclusions. The electrostatic component of that exception is
# already handled (since a scaling factor *can* represent the full
# flexibility of electrostatic exceptions). The vdW component must be
# handled as 1-4 A- and B-coefficients. If vdW types are too compressed
# for this to be handled correctly, use "addLJType" to expand the types
# by 1 (so the tables only get as big as they *need* to get, to make
# sure they continue to fit in CUDA shared memory).
#
# The way we do this is to fill all of the elements with None, then fill
# them in as we walk through the 1-4 exceptions. If we hit a case where
# the target element is not None and *doesn't* equal the computed A- and
# B-coefficients, we have to expand our type list (assume all type
# names will have the same L-J parameters, which has been a fair
# assumption in my experience).
if not self.adjusts: return
for i in range(len(data['LENNARD_JONES_14_ACOEF'])):
data['LENNARD_JONES_14_ACOEF'][i] = None
data['LENNARD_JONES_14_BCOEF'][i] = None
ii = 0
replaced_atoms = set()
while True:
needed_split = False
for pair in self.adjusts:
a1, a2 = pair.atom1, pair.atom2
i, j = sorted([a1.nb_idx - 1, a2.nb_idx - 1])
idx = data['NONBONDED_PARM_INDEX'][ntypes*i+j] - 1
eps = pair.type.epsilon
rmin = pair.type.rmin
rmin6 = rmin * rmin * rmin * rmin * rmin * rmin
acoef = eps * rmin6*rmin6
bcoef = 2 * eps * rmin6
if data['LENNARD_JONES_14_ACOEF'][idx] is not None:
if abs(data['LENNARD_JONES_14_ACOEF'][idx] - acoef) > SMALL:
# Need to split out another type
needed_split = True
assert a1 not in replaced_atoms or a2 not in replaced_atoms
# Only add each atom as a new type ONCE
if a1 in replaced_atoms:
mask = '@%d' % (a2.idx+1)
replaced_atoms.add(a2)
else:
mask = '@%d' % (a1.idx+1)
replaced_atoms.add(a1)
addLJType(self, mask, radius_14=0, epsilon_14=0).execute()
ntypes += 1
# None-out all of the added terms
j = ntypes - 1
for i in range(j):
idx2 = data['NONBONDED_PARM_INDEX'][ntypes*i+j] - 1
data['LENNARD_JONES_14_ACOEF'][idx2] = None
data['LENNARD_JONES_14_BCOEF'][idx2] = None
# We can stop here, since the next loop through the
# explicit exclusions will fill this in
else:
data['LENNARD_JONES_14_ACOEF'][idx] = acoef
data['LENNARD_JONES_14_BCOEF'][idx] = bcoef
ii += 1
if not needed_split:
break
# The following should never happen
assert ii <= len(self.atoms)+1, 'Could not resolve all exceptions. ' \
'Some unexpected problem with the algorithm'
# Now go through and change all None's to 0s, as these terms won't be
# used for any exceptions, anyway
for i, item in enumerate(data['LENNARD_JONES_14_ACOEF']):
if item is None:
assert data['LENNARD_JONES_14_BCOEF'][i] is None, \
'A- and B- coefficients must be in lock-step!'
data['LENNARD_JONES_14_ACOEF'][i] = 0.0
data['LENNARD_JONES_14_BCOEF'][i] = 0.0
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[docs]def ConvertFromPSF(struct, params, title=''):
"""
This function instantiates a ChamberParm instance from a data structure
instantiated by a CHARMM PSF.
Parameters
----------
struct : Structure
The structure object (typically loaded from a PSF file)
params : CharmmParameterSet
The parameter set describing the parameters of the input struct
title : str=''
The title to assign to the topology file
Returns
-------
ChamberParm
ChamberParm instance with all parameters loaded
"""
# Make sure all atom types are strings, not integers
int_starting = str(struct.atoms[0].atom_type) != struct.atoms[0].type
if int_starting:
for atom in struct.atoms:
atom.type = str(atom.atom_type)
parm = ChamberParm.from_structure(struct)
parm.parm_data['FORCE_FIELD_TYPE'] = fftype = []
if params.parametersets == []:
params.parametersets.append('')
for pset in params.parametersets:
if 'CHARMM' not in pset: # needed to trigger "charmm_active"...
pset = 'CHARMM: %s' % pset
fftype.extend([len(params.parametersets), pset])
# Convert atom types back to integers if that's how they started
if int_starting:
for atom in struct.atoms:
atom.type = int(atom.atom_type)
return parm