Skip to content
Snippets Groups Projects
Commit a4162b1d authored by cmaffeo2's avatar cmaffeo2
Browse files

Added ENM for atomistic structures

parent c5253f62
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
atomicModel.py
+ 243
− 3
View file @ a4162b1d
......@@ -200,6 +200,7 @@ class atomicModel():
self.nodes = None
self.nodeTypeCounts = None
self.ntScale = ntScale
self.part = part
# self._nbParams = set()
self._bondParams = set()
......@@ -214,6 +215,7 @@ class atomicModel():
random.seed(randomSeed)
self._buildModel(part)
self.latticeType = self._determineLatticeType(part)
def _helixStrandsToEnds(self, helixStrands):
......@@ -247,7 +249,7 @@ class atomicModel():
origins = []
if props.get('point_type') == PointType.ARBITRARY:
# TODO add code to encode Parts with ARBITRARY point configurations
raise NotImplementedError("Not implemented")
raise Exception("Not implemented")
else:
vh_props, origins = part.helixPropertiesAndOrigins()
# print(' VIRTUAL HELICES:', vh_props)
......@@ -436,6 +438,29 @@ class atomicModel():
node1.addXover(node2, isFwd1)
node2.addXover(node1, isFwd2)
def _determineLatticeType(self, part):
props, origins = part.helixPropertiesAndOrigins()
origins = [np.array(o) for o in origins]
neighborVecs = []
for i in range(len(props)):
js = [int(j.strip('[,]')) for j in props['neighbors'][i].split()]
for j in js:
neighborVecs.append( origins[j]-origins[i] )
dots = [nv2.dot(nv1)/(np.linalg.norm(nv1)*np.linalg.norm(nv2)) for nv1 in neighborVecs for nv2 in neighborVecs]
nDots = len(dots)
nAngled = np.sum( np.abs(dots)-0.5 < 0.2 )
if float(nAngled)/float(nDots) == 0:
return 'square'
elif float(nAngled)/float(nDots) < 0.05:
print('WARNING: unusual neighbors in square lattice structure')
return 'square'
elif float(nAngled)/float(nDots) > 0.4:
return 'honeycomb'
else:
raise Exception("Could not identify lattice type")
def backmap(self, simplerModel, simplerModelCoords,
dsDnaHelixNeighborDist=50, dsDnaAllNeighborDist=30,
......@@ -639,6 +664,11 @@ class atomicModel():
# -------------------------- #
# Methods for printing model #
# -------------------------- #
def writeNamdFiles(self,prefix):
self.writePdbPsf(prefix)
self.writeENM(prefix)
self.writeNamdFile(prefix)
def writePdbPsf(self, prefix):
bonds,angles,dihedrals,impropers = [[],[],[],[]]
with open("%s.pdb" % prefix,'w') as pdb, \
......@@ -773,5 +803,215 @@ class atomicModel():
psf.write("\n\n 1 0 !NGRP\n\n")
def writeENM(self):
return
def writeENM(self, prefix):
if self.latticeType == "square":
enmTemplate = enmTemplateSQ
elif self.latticeType == "honeycomb":
enmTemplate = enmTemplateHC
else:
raise Exception("Lattice type '%s' not supported" % self.latticeType)
with open("%s.exb" % prefix,'w') as fh:
# natoms=0
for seg in self.segments:
for nt1 in seg:
other = []
nt2 = nt1.basepair
if nt2 is not None:
if nt2.firstAtomicIndex > nt1.firstAtomicIndex:
other.append((nt2,'pair'))
nt2 = nt2.ntAt3prime
if nt2 is not None:
other.append((nt2,'paircross'))
nt2 = nt1.ntAt3prime
if nt2 is not None:
other.append((nt2,'stack'))
nt2 = nt2.basepair
if nt2 is not None:
other.append((nt2,'cross'))
a1 = nt1.atoms(transform=False)
for nt2,key in other:
key = ','.join((key,nt1.sequence[0],nt2.sequence[0]))
for n1, n2, d in enmTemplate[key]:
a2 = nt2.atoms(transform=False)
try:
i,j = [a.index(n) for a,n in zip((a1,a2),(n1,n2))]
except:
continue
k = 0.1
if self.latticeType == 'honeycomb':
correctionKey = ','.join(key,n1,n2)
dk,dr = enmCorrectionsHC[correctionKey]
k += dk
d += dr
i += nt1.firstAtomicIndex
j += nt2.firstAtomicIndex
fh.write("bond %d %d %.2f %.2f\n".format(i,j,k,d))
props, origins = self.part.helixPropertiesAndOrigins()
origins = [np.array(o) for o in origins]
## Push bonds
ntDict = dict()
xo = dict()
helixIds = [hid for s in self.segments for hid in s.nodes.keys()]
helixIds = sorted(list(set(helixIds)))
for h1 in helixIds:
xo[h1] = dict()
for h2 in helixIds:
xo[h1][h2] = []
for seg in self.segments:
for nt in seg:
hid = nt.prop['helixId']
# if hid not in ntDict:
# ntDict[hid] = dict()
if hid not in xo:
xo[hid] = dict()
## Add nt to ntDict
zidx = nt.prop['zIdx']
isFwd = nt.prop['isFwd']
ntDict[(hid,zidx,isFwd)] = nt
## Find crossovers
for nt2 in (nt.ntAt5prime, nt.ntAt3prime):
if nt2 is not None:
hid2 = nt2.prop['helixId']
if hid2 != hid:
xo[hid][hid2].append(zidx)
props = self.part.getModelProperties()
if props.get('point_type') == PointType.ARBITRARY:
raise Exception("Not implemented")
else:
props, origins = self.part.helixPropertiesAndOrigins()
neighborHelices = dict()
for hid in helixIds:
neighborHelices[hid] = [int(j.strip('[,]')) for j in props['neighbors'][hid].split()]
if self.latticeType == 'honeycomb':
minDist = 7
elif self.latticeType == 'square':
minDist = 10
for seg in self.segments:
for nt in seg:
hid = nt.prop['helixId']
zidx = nt.prop['zIdx']
isFwd = nt.prop['isFwd']
for hid2 in neighborHelices[hid]:
dist = min([z-zidx for z in xo[hid][hid2]])
key = (hid2,zidx,isFwd)
if dist >= minDist and key in ntDict:
nt2 = ntDict[key]
i = nt.firstAtomicIndex
j = nt2.firstAtomicIndex
if j > i:
i += nt.atoms(transform=False).index('P')
j += nt2.atoms(transform=False).index('P')
fh.write("bond %d %d %.2f %.2f\n" % (i,j,1,31))
def writeNamdFile(self,prefix):
with open("%s.namd" % prefix,'w') as fh:
fh.write("""#############################################################
## CONFIGURATION FILE FOR ORIGAMI STRUCTURE PREDICTION ##
#############################################################
#############################################################
## ADJUSTABLE PARAMETERS ##
#############################################################\n""")
fh.write("set prefix %s\n" % prefix)
fh.write("""set nLast 0; # increment when continueing a simulation
set n [expr $nLast+1]
set out output/$prefix-$n
set temperature 300
structure $prefix.psf
coordinates $prefix.pdb
outputName $out
XSTfile $out.xst
DCDfile $out.dcd
#############################################################
## SIMULATION PARAMETERS ##
#############################################################
# Input
paraTypeCharmm on
parameters charmm36.nbfix/par_all36_na.prm
parameters charmm36.nbfix/par_water_ions_na.prm
wrapAll off
# Force-Field Parameters
exclude scaled1-4
1-4scaling 1.0
switching on
switchdist 8
cutoff 10
pairlistdist 12
margin 30
# Integrator Parameters
timestep 2.0 ;# 2fs/step
rigidBonds all ;# needed for 2fs steps
nonbondedFreq 1
fullElectFrequency 3
stepspercycle 12
# PME (for full-system periodic electrostatics)
PME no
PMEGridSpacing 1.2
# Constant Temperature Control
langevin on ;# do langevin dynamics
# langevinDamping 1 ;# damping coefficient (gamma); used in original study
langevinDamping 0.1 ;# less friction for faster relaxation
langevinTemp $temperature
langevinHydrogen off ;# don't couple langevin bath to hydrogens
# output
useGroupPressure yes
xstFreq 4800
outputEnergies 4800
dcdfreq 4800
restartfreq 48000
#############################################################
## EXTRA FORCES ##
#############################################################
# ENM and intrahelical extrabonds
extraBonds on
extraBondsFile $prefix.exb
#############################################################
## RUN ##
#############################################################
# Continuing a job from the restart files
cellBasisVector1 1000 0 0
cellBasisVector2 0 1000 0
cellBasisVector3 0 0 1000
if {$nLast == 0} {
temperature 300
minimize 4800
} else {
bincoordinates output/$prefix-$nLast.restart.coor
binvelocities output/$prefix-$nLast.restart.vel
}
run 96000000
""")
charmm36.nbfix/par_all36_na.prm 0 → 100644
+ 1253
− 0
View file @ a4162b1d
This diff is collapsed.
charmm36.nbfix/par_water_ions_na.prm 0 → 100644
+ 203
− 0
View file @ a4162b1d
!read para card flex @app
!* Parameters for water and ions
!*
!ATOMS
!MASS 1 HT 1.00800 ! TIPS3P WATER HYDROGEN
!MASS 2 HX 1.00800 ! hydroxide hydrogen
!MASS 3 OT 15.99940 ! TIPS3P WATER OXYGEN
!MASS 4 OX 15.99940 ! hydroxide oxygen
!MASS 5 LIT 6.94100 ! Lithium ion
!MASS 6 SOD 22.98977 ! Sodium Ion
!MASS 7 MG 24.30500 ! Magnesium Ion
!MASS 8 POT 39.09830 ! Potassium Ion
!MASS 9 CAL 40.08000 ! Calcium Ion
!MASS 10 RUB 85.46780 ! Rubidium Ion
!MASS 11 CES 132.90545 ! Cesium Ion
!MASS 12 BAR 137.32700 ! Barium Ion
!MASS 13 ZN 65.37000 ! zinc (II) cation
!MASS 14 CAD 112.41100 ! cadmium (II) cation
!MASS 15 CLA 35.45000 ! Chloride Ion
!MASS 16 OTMG 15.99940 ! TIPS3P WATER OXYGEN by jejoong
BONDS
!
!V(bond) = Kb(b - b0)**2
!
!Kb: kcal/mole/A**2
!b0: A
!
!atom type Kb b0
!
HT HT 0.0 1.5139 ! from TIPS3P geometry (for SHAKE w/PARAM)
HT OT 450.0 0.9572 ! from TIPS3P geometry
OX HX 545.0 0.9700 ! hydroxide ion
!jejoong
MG OTMG 10000.00 1.9400 !
OTMG HT 450.000 0.9572
ANGLES
!
!V(angle) = Ktheta(Theta - Theta0)**2
!
!V(Urey-Bradley) = Kub(S - S0)**2
!
!Ktheta: kcal/mole/rad**2
!Theta0: degrees
!Kub: kcal/mole/A**2 (Urey-Bradley)
!S0: A
!
!atom types Ktheta Theta0 Kub S0
!
HT OT HT 55.0 104.52 ! FROM TIPS3P GEOMETRY
!jejoong
HT OTMG HT 55.000 104.5200
HT OTMG MG 0.000 104.5200
OTMG MG OTMG 1000.000 90.0000
DIHEDRALS
!
!V(dihedral) = Kchi(1 + cos(n(chi) - delta))
!
!Kchi: kcal/mole
!n: multiplicity
!delta: degrees
!
!atom types Kchi n delta
!
HT OTMG MG OTMG 0.000 6 180.0000
!
IMPROPER
!
!V(improper) = Kpsi(psi - psi0)**2
!
!Kpsi: kcal/mole/rad**2
!psi0: degrees
!note that the second column of numbers (0) is ignored
!
!atom types Kpsi psi0
!
NONBONDED nbxmod 5 atom cdiel shift vatom vdistance vswitch -
cutnb 14.0 ctofnb 12.0 ctonnb 10.0 eps 1.0 e14fac 1.0 wmin 1.5
!TIP3P LJ parameters
HT 0.0 -0.046 0.2245
OT 0.0 -0.1521 1.7682
!jejoong
OTMG 0.000000 -0.152100 1.768200 ! TIP3P oxygen
!for hydroxide
OX 0.000000 -0.120000 1.700000 ! ALLOW POL ION
! JG 8/27/89
HX 0.000000 -0.046000 0.224500 ! ALLOW PEP POL SUL ARO ALC
! same as TIP3P hydrogen, adm jr., 7/20/89
!ions
LIT 0.0 -0.00233 1.2975 ! Lithium
! From S Noskov, target ddG(Li-Na) was 23-26.0 kcal/mol (see JPC B, Lamoureux&Roux,2006)
SOD 0.0 -0.0469 1.41075 ! new CHARMM Sodium
! ddG of -18.6 kcal/mol with K+ from S. Noskov
MG 0.0 -0.0150 1.18500 ! Magnesium
! B. Roux dA = -441.65
POT 0.0 -0.0870 1.76375 ! Potassium
! D. Beglovd and B. Roux, dA=-82.36+2.8 = -79.56 kca/mol
CAL 0.0 -0.120 1.367 ! Calcium
! S. Marchand and B. Roux, dA = -384.8 kcal/mol
RUB 0.0000 -0.15 1.90 ! Rubidium
! delta A with respect to POT is +6.0 kcal/mol in bulk water
CES 0.0 -0.1900 2.100 ! Cesium
! delta A with respect to POT is +12.0 kcal/mol
BAR 0.0 -0.150 1.890 ! Barium
! B. Roux, dA = dA[calcium] + 64.2 kcal/mol
ZN 0.000000 -0.250000 1.090000 ! Zinc
! RHS March 18, 1990
CAD 0.000000 -0.120000 1.357000 ! Cadmium
! S. Marchand and B. Roux, from delta delta G
CLA 0.0 -0.150 2.27 ! Chloride
! D. Beglovd and B. Roux, dA=-83.87+4.46 = -79.40 kcal/mol
NBFIX
! Emin Rmin
! (kcal/mol) (A)
!SOD CLA -0.083875 3.731 ! From osmotic pressure calibration, J. Phys.Chem.Lett. 1:183-189
!POT CLA -0.114236 4.081 ! From osmotic pressure calibration, J. Phys.Chem.Lett. 1:183-189
! monovalent cation -- protein acetate (OC), lipid acetate (OCL)
!LIT OC -0.016721 2.9975 ! standard
!LIT OC -0.016721 3.1875 ! +0.19
!LIT OCL -0.016721 3.1875 ! +0.19
!SOD OC -0.075020 3.11075 ! standard
!SOD OC -0.075020 3.20075 ! +0.09
!SOD OCL -0.075020 3.20075 ! +0.09
!POT OC -0.102176 3.46375 ! standard
!POT OC -0.102176 3.54375 ! +0.08
!POT OCL -0.102176 3.54375 ! +0.08
! monovalent cation -- lipid phosphate (O2L), nucleic acid phosphate (ON3)
!LIT OC -0.016721 2.9975 ! standard
!LIT O2L -0.016721 3.1875 ! +0.19
LIT ON3 -0.016721 3.1875 ! +0.19
!SOD OC -0.075020 3.11075 ! standard
!SOD O2L -0.075020 3.20075 ! +0.09
SOD ON3 -0.075020 3.20075 ! +0.09
!POT OC -0.102176 3.46375 ! standard
!POT O2L -0.102176 3.54375 ! +0.08
POT ON3 -0.102176 3.54375 ! +0.08
! monovalent cation -- chloride
!LIT CLA -0.018695 3.5675 ! standard
LIT CLA -0.018695 3.7975 ! +0.23
!SOD CLA -0.083875 3.68075 ! standard
SOD CLA -0.083875 3.74075 ! +0.06
!POT CLA -0.114237 4.03375 ! standard
POT CLA -0.114237 4.08375 ! +0.05
! magnesium -- chloride
!OTMG CLA -0.151046 4.03820 ! standard
OTMG CLA -0.151046 4.05320 ! +0.015
! magnesium -- acetate / phosphate
!OTMG OC -0.135099 3.46820 ! standard
!OTMG OC -0.135099 3.52320 ! +0.055
!OTMG OCL -0.135099 3.52320 ! +0.055
!OTMG O2L -0.135099 3.52320 ! +0.055
OTMG ON3 -0.135099 3.52320 ! +0.055
END
! The following section contains NBFixes for sodium interacting with
! carboxylate oxygens of various CHARMM force fields. It will generate
! level -1 warnings whenever any of these force fields have not been
! read prior to the current stream file. Since we don't want to force
! the user to always read all the force fields, we're suppressing the
! warnings. The only side effect is that you will have "most severe
! warning was at level 0" at the end of your output. Also note that
! the user is responsible for reading the current file last if they
! want the NBFixes to apply. A more elegant solution would require new
! features to be added to CHARMM.
! parallel fix, to avoid duplicated messages in the log
!set para
!if ?NUMNODE gt 1 set para node 0
!set wrn ?WRNLEV
! Some versions of CHARMM don't seem to initialize wrnlev...
!if "@WRN" eq "?WRNLEV" set wrn 5
!set bom ?bomlev
!WRNLEV -1 @PARA
!BOMLEV -1 @PARA
!!read para card flex append
!* NBFix between carboxylate and sodium
!*
! These NBFixes will only apply if the main files have been read in first!!!
!NBFIX
!SOD OC -0.075020 3.190 ! For prot carboxylate groups
!SOD OCL -0.075020 3.190 ! For lipid carboxylate groups
!SOD OC2D2 -0.075020 3.190 ! For carb carboxylate groups
!SOD OG2D2 -0.075020 3.190 ! For CGenFF carboxylate groups
!END
!BOMLEV @bom @PARA
!WRNLEV @wrn @PARA
return
charmm36.nbfix/top_all36_na.rtf 0 → 100644
+ 853
− 0
View file @ a4162b1d
This diff is collapsed.
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment