Use CHARMM force field for PEG.

This is part of my recent research. To make it simple, I am going to build a TIP3P water box with one PEG polymer (5 units).

 

CHARMM force field for ethers has been well testified in a the following paper with experimental data. So I have to say it the best force field for PEG by now.

 

Lee, Hwankyu, Richard M. Venable, Alexander D. MacKerell, and Richard W. Pastor. “Molecular Dynamics Studies of Polyethylene Oxide and Polyethylene Glycol: Hydrodynamic Radius and Shape Anisotropy.” Biophysical Journal 95, no. 4 (August 15, 2008): 1590–99. doi:10.1529/biophysj.108.133025.
 ————————————————-
You’d better know something about MD simulation or CHARMM. What you need to do is:
1. Download force field files.
2. Construct PEG in vacuum.
3. Build a water box.
Is it easy? This cost my four days to figure out how to do that. Dr. Justin Lemkul gave me a lot of help.
———————————-
The first part, force field files.
You need to download the two files (RTF, PRM) at https://www.charmm.org/ubbthreads/ubbthreads.php?ubb=showflat&Number=7023&page=2.
It is unbelievable that the published version for ethers missed some important termini for PEG. And you need to modify the par_all35_ethers-oh.prm file, because there are some duplicate terms. I choose to comment lines 157 to 166, which are old parameters. If you don’t do that. You will find that CHARMM will give error about these three dihedrals.
———————————-
The second part, build PEG.
I knew nothing about CHARMM scripts. So I tried to figure out how CHARMMing.org and CHARMM-gui work. These online tools are really awesome! And I used the scripts got from them to do this step. The script is below.
* make a PDB file out of a sequence
*

read rtf card name top_all35_ethers-oh.rtf
read param card name par_all35_ethers-oh.prm

read sequ card

4
PEGM PEGM PEGM PEGM

gene a setup first GCL0 last GCL3

ic param
ic seed 1 C1 1 O1 1 C2
ic build

hbuild
coor orient

write coor pdb name peg5vac.pdb
* PDB file of custom sequence.
*
write psf card name peg5vac.psf
* PSF
*

stop
—————————–
The third part, water box.
At the beginning, I tried to use scripts from CHARMMing and CHARMM-gui, but even I uploaded the ethers rtf and prm files, there would be errors. After trying hundreds of modification of these scripts, I gave it up.
Unintentionally, I found VMD can do that in a very easy way.
I am using iMac. So load the files in VMD by following commands using files got from the second part.
vmd peg5vac.psf peg5vac.pdb
Open tcl/tk console and input the following commands. New psf and pdb files will be generated. The size and center of box can be used in the simulation.
package require solvate
solvate peg5vac.psf peg5vac.pdb -t 15 -o peg5wat
set everyone [atomselect top all]
measure minmax $everyone
measure center $everyone
 Commands are cited from http://www.ks.uiuc.edu/Training/Tutorials/namd/namd-tutorial-unix-html/node7.html.
————————————–
Now you have peg5wat.psf and peg5wat.pdb. You can run your simulation!

Truncated octahedral cell in NAMD.

I used cubic box before, but it had too many water molecules for big system. Truncated octahedral could reduce 50% size of the overall system.

I got two solutions for this. To visualize the box correctly, add wrapall on and wrapnearest on in your configuration file.

—————————————————–

The first solution was got from http://nefeli.mbg.duth.gr/cgi-bin/wiki.pl?action=browse;diff=1;id=Performing_a_molecular_dynamics_simulation_in_a_truncated_octahedral_cell.

In the NAMD scripts you will have to define the unit cell translations for the truncated octahedron (the cellbasis keywords). The matrix you need is

d 0 0
0 d 0
d/2 d/2 d/2

where d is the length of the edge of the cube you’ve specified in the vmd script. For the example shown above you would have to specify

cellBasisVector1        32.00    0.00    0.00   
cellBasisVector2         0.00   32.00    0.00   
cellBasisVector3        16.00   16.00   16.00   
cellOrigin               0.00    0.00    0.00

——————————————————–

The second solution was got form http://ambermd.org/namd/namd_amber.html.

Truncated Octahedron (TO)

To use a TO in sander/pmemd you have to use the solvateOct command in tleap/sleapm, or to modify the last line of the crd/rst file to the corresponding values. NAMD will not read the box information from the crd/rst file generated from tleap/sleap. Instead you will have to specify the box information using the cellBasisVector1, 2, 3. If the d is the length of the edge of the TO, than the three vectors should be:
cellBasisVector1  d         0.0            0.0
cellBasisVector2  (-1/3)*d (2/3)sqrt(2)*d  0.0
cellBasisVector3  (-1/3)*d (-1/3)sqrt(2)*d (-1/3)sqrt(6)*d
The last line of your crd/rst file gnerated from tleap/sleap or from sander/pmemd should look like:
d d d 109.4712190 109.4712190 109.4712190
------------------------------------------------------

Fix atoms in simulation with NAMD

I am working on NAMD recently. I used AMBER force field files. NAMD did have much better parallel computing for big system.

The solution was got from http://www.ks.uiuc.edu/Research/namd/mailing_list/namd-l.2009-2010/3969.html.

1) you need 1st create a new .pdb file with B col specified as ‘0’ for unaffected atoms, and ‘1’ for fixed or constrained atoms. you can do it with VMD.
after you load your protein and open the TK console, type the following for example,
set all [atomselect top all]
set fixatom [atomselect top “protein”]
$all set beta 0
$fixatom set beta 1
$all writepdb fixprotein.pdb

2) then you simply implement the fix or constrain patch into your .conf file
a short example,
fixedAtoms on
fixedAtomsFile fixprotein.pdb
fixedAtomsCol B