MBG wiki | RecentChanges | Blog | 2021-09-24 | 2021-09-23

Running a molecular dynamics simulation of a tetrapeptide in a truncated octahedral cell



Go read & do the previous tutorials on Rop and protein G, then come back … (read also the page on truncated octahedra)



If you are seriously into peptides, do read Spoel's thesis



Peptide PDB generation

The program you need is ribosome. Prepare a script file with your sequence :

# cat > ribosome.script
title FPGF
default helix
res phe
res pro
res gly
res phe
<CTRL-D>

Run the program :

/usr/local/bin/ribosome ribosome.script FPGF.pdb /usr/local/bin/ribosome.dat

This should create a PDB file with the name FPGF.pdb containing your tetrapeptide. The backbone conformation is helical (as requested in the script file), the side chains' conformations are random.

Continuing ...

Use moleman to align the axes of inertia :

#!/bin/tcsh -f

#
# This will read a PDB file and rotate/translate it so that
#  1. the centre of gravity will be at the origin
#  2. the axes of inertia will be aligned with the orthogonal frame
#

lx_moleman2 >& moleman.log << eof
/usr/local/xutil/moleman2.lib
REad FPGF.pdb
XYz ALign_inertia_axes
WRite aligned.pdb
quit
eof

exit

Edit the PDB file created by moleman, remove the header, give all atoms a chain identifier (e.g. 'A').

Use psfgen to generate PSF and PDB files :

#!/bin/tcsh -f

/usr/local/NAMD_2.5/psfgen >& psfgen.log << END
topology top_all27_prot_na.inp

segment A {
 pdb aligned.pdb
 }
alias atom ILE CD1 CD
coordpdb aligned.pdb A


guesscoord
writepsf psfgen.psf
writepdb psfgen.pdb

END

exit

Use pdbset to determine limits (along x, y and z) of the aligned molecule. Decide the size of the cube of water that has to be added.

Use vmd to add waters and ions for a truncated octahedral cell. Your vmd script should look similar to

#
# Make water box
#
package require vexpr
package require toctsolvate
toctsolvate psfgen.psf psfgen.pdb -o hydrated -minmax {{-17 -17 -17} {17 17 17}}  -b 1.80 

#
# Add ions to neutralise charge
#
package require autoionize
autoionize -psf hydrated.psf -pdb hydrated.pdb -is 0.20

#
# Prepare restraints files
#
mol load psf ionized.psf pdb ionized.pdb

set all [atomselect top all]
set sel [atomselect top "protein and name CA"]
$all set beta 0
$sel set beta 0.5
$all writepdb restrain_ca.pdb


set all [atomselect top all]
set to_fix [atomselect top "protein and backbone"]
$all set beta 0
$to_fix set beta 1
$all writepdb fix_backbone.pdb

Examine the system (ionised.pdb).

Heat it up (for a real run you should increase the time spend in the heating phase) :

#
# Input files
#
structure               ionized.psf
coordinates             ionized.pdb
parameters              par_all27_prot_na.inp
paraTypeCharmm          on

#
# Output files & writing frequency for DCD
# and restart files
#
outputname              output/heat_out
binaryoutput            off
restartname             output/restart
restartfreq             1000
binaryrestart           yes
dcdFile                 output/heat_out.dcd
dcdFreq                 200
DCDunitcell             on

#
# Frequencies for logs and the xst file
#
outputEnergies          20
outputTiming            200
xstFreq                 200

#
# Timestep & friends
#
timestep                2.0
stepsPerCycle           8
nonBondedFreq           2
fullElectFrequency      4

#
# Simulation space partitioning
#
switching               on
switchDist              8
cutoff                  9
pairlistdist            11

#
# Basic dynamics
#
temperature             0
COMmotion               no
dielectric              1.0
exclude                 scaled1-4
1-4scaling              1.0
rigidbonds              all

#
# Particle Mesh Ewald parameters. 
#
Pme                     on
PmeGridsizeX            32                      # <===== CHANGE ME
PmeGridsizeY            32                      # <===== CHANGE ME
PmeGridsizeZ            27                      # <===== CHANGE ME


#
# Periodic boundary things
#
wrapWater               on
wrapNearest             on

cellBasisVector1        32.00    0.00    0.00   # <===== CHANGE ME
cellBasisVector2         0.00   32.00    0.00   # <===== CHANGE ME
cellBasisVector3        16.00   16.00   16.00   # <===== CHANGE ME
cellOrigin               0.00    0.00    0.00   # <===== CHANGE ME

#
# Fixed atoms for initial heating-up steps
#
fixedAtoms              on
fixedAtomsForces        on
fixedAtomsFile          fix_backbone.pdb
fixedAtomsCol           B

#
# Restrained atoms for initial heating-up steps
#
constraints             on
consRef                 restrain_ca.pdb
consKFile               restrain_ca.pdb
consKCol                B

#
# Langevin dynamics parameters
#
langevin                on
langevinDamping         10
langevinTemp            320                     # <===== Check me
langevinHydrogen        on

langevinPiston          on
langevinPistonTarget    1.01325
langevinPistonPeriod    200
langevinPistonDecay     100
langevinPistonTemp      320                     # <===== Check me

useGroupPressure        yes


##########################################
# The actual minimisation and heating-up 
# protocol follows. The number of steps
# shown below are too small for a real run
##########################################

#
# run one step to get into scripting mode
#
minimize                0

#
# turn off pressure control until later
#
langevinPiston          off

#
# minimize nonbackbone atoms
#
minimize                400                     ;# <===== CHANGE ME
output                  output/min_fix

#
# min all atoms
#
fixedAtoms              off
minimize                400                     ;# <===== CHANGE ME
output                  output/min_all

#
# heat with CAs restrained
#
set temp 20;
while { $temp < 321 } {                         ;# <===== Check me
langevinTemp            $temp
run                     400                     ;# <===== CHANGE ME
output                  output/heat_ca
set temp [expr $temp + 20]                      
}

#
# equilibrate volume with CAs restrained
#
langevinPiston          on
run                     400                     ;# <===== CHANGE ME
output                  output/equil_ca

#
# equilibrate volume without restraints
#
constraintScaling       0
run                     2000                    ;# <===== CHANGE ME

Production run :

#
# Input files
#
structure               ionized.psf
coordinates             heat_out.coor
velocities              heat_out.vel
extendedSystem          heat_out.xsc
parameters              par_all27_prot_na.inp
paraTypeCharmm          on

#
# Output files & writing frequency for DCD
# and restart files
#
outputname              output/equi_out
binaryoutput            off
restartname             output/restart
restartfreq             1000
binaryrestart           yes
dcdFile                 output/equi_out.dcd
dcdFreq                 200
DCDunitcell             on

#
# Frequencies for logs and the xst file
#
outputEnergies          20
outputTiming            200
xstFreq                 200

#
# Timestep & friends
#
timestep                2.0
stepsPerCycle           8
nonBondedFreq           2
fullElectFrequency      4

#
# Simulation space partitioning
#
switching               on
switchDist              8
cutoff                  9
pairlistdist            11

#
# Basic dynamics
#
COMmotion               no
dielectric              1.0
exclude                 scaled1-4
1-4scaling              1.0
rigidbonds              all

#
# Particle Mesh Ewald parameters. 
#
Pme                     on
PmeGridsizeX            32                      # <===== CHANGE ME
PmeGridsizeY            32                      # <===== CHANGE ME
PmeGridsizeZ            27                      # <===== CHANGE ME


#
# Periodic boundary things
#
wrapWater               on
wrapNearest             on
wrapAll                 on


#
# Langevin dynamics parameters
#
langevin                on
langevinDamping         1
langevinTemp            320                     # <===== Check me
langevinHydrogen        on

langevinPiston          on
langevinPistonTarget    1.01325
langevinPistonPeriod    200
langevinPistonDecay     100
langevinPistonTemp      320                     # <===== Check me

useGroupPressure        yes

firsttimestep           9600                    # <===== CHANGE ME
run                     5000000                ;# <===== CHANGE ME