* In this file we will build an alanine hexa-peptide * in an a-helical conformation and perform 10 ps * of molecular dynamics on the molecule. This * will illustrate the basic features of setting-up * molecular dynamics calculations as well as * building polypeptides in specific conformations. * ! The first task is to read in the topology and parameters for ! the systems we want to simulate. This we will accomplish in ! manner similar to that done in the previous input. ! First we'll read in the topology file, this is the same as last ! time except we'll change the IC's so the structure is built in ! an a-helical conformation read rtf card ! translated - read Residue Topology File using an * TOPOLOGY FILE FOR PROTEINS * USING EXPLICIT HYDROGEN ATOMS: VERSION 19 * 20 1 ! Version number ! atm #, atm name, atm mass MASS 1 H 1.00800 ! hydrogen which can h-bond to neutral atom MASS 2 HC 1.00800 ! - " - charged atom MASS 11 C 12.01100 ! carbonyl carbon MASS 12 CH1E 13.01900 ! extended atom carbon w/ one hydrogen MASS 14 CH3E 15.03500 ! - " - three MASS 38 NH1 14.00670 ! peptide nitrogen bound to one hydrogen MASS 40 NH3 14.00670 ! nitrogen bound to three hydrogens MASS 51 O 15.99940 ! carbonyl oxygen MASS 52 OC 15.99940 ! carboxy oxygen DECL -C ! declarations are useful to describe connectivity to DECL -O ! residues before and after in the chain DECL +N DECL +H DECL +CA AUTOGENERATE ANGLES ! can specify autogenerate angles to make life simple DEFAult_patches FIRSt NTERminal LAST CTERminal ! can also specify caps ! We will now define the alanine residue by a residue and two patches, ! by doing it in this way we are allowed to build-up multiple residue ! chains which are then "capped" with the end patches. ! The alanine we will build is: ! ! CB O ! | || ! H-N--CA--C ! ! We will then patch on the other H and O RESI ALA 0.00000 ! name the residue and give the overall charge GROU ATOM N NH1 -0.35 ! specify all atoms, you can group them in ATOM H H 0.25 ! accordance to neutral groups as done here ATOM CA CH1E 0.10 ! the atom specification consists of the name ! of the atom in the specific residue, its type and its atomic partial ! charge in the specific group GROU ATOM CB CH3E 0.00 ! the alanine sidechain GROU ATOM C C 0.55 ! now the carbonyl group ATOM O O -0.55 ! Next we specify all bonds making up the bonded topology BOND N CA CA C C +N C O N H BOND CA CB ! And finally the dihedrals and impropers, one dihedral is defined ! around every rotable bond and impropers are defined to maintain ! chirality when explicit C-H hydrogens are not present, or sometimes ! to keep rings planer. DIHE -C N CA C N CA C +N CA C +N +CA IMPH N -C CA H C CA +N O CA N C CB ! Lastly we can specify a given set of internal coordinates (ICs) which ! permit us to build the structure without reference to an external set ! of atomic coordinates. The ICs specify bonds, angles and dihedrals ! (both proper and improper), however, we'll just give the desired ! dihedrals and "fill" the table for bonds and angles from the parameters ! we'll input below. IC -C CA *N H 0.0000 0.00 180.00 0.00 0.0000 !*** We will set phi, psi to -60, -60 to be near the right-handed helix IC -C N CA C 0.0000 0.00 -60.00 0.00 0.0000 ! phi IC N CA C +N 0.0000 0.00 -60.00 0.00 0.0000 ! psi !******************************************************************* IC +N CA *C O 0.0000 0.00 180.00 0.00 0.0000 IC CA C +N +CA 0.0000 0.00 180.00 0.00 0.0000 IC N C *CA CB 0.0000 0.00 120.00 0.00 0.0000 ! We now include topologies for the "patch residues", first the ! N-terminal patch. Again, this is specified much like the ! regular residue, except we can do some wierd things like delete ! atoms and replace them by others. This is done for the amide ! hydrogen, whose "name" (H) no longer exists in the zwitter ion form PRES NTER 1.00000 ! total charge of patch residue nter is 1.0 GROU ATOM HT1 HC 0.35 ATOM HT2 HC 0.35 ATOM N NH3 -0.30 ATOM HT3 HC 0.35 DELETE ATOM H ATOM CA CH1E 0.25 BOND HT1 N HT2 N HT3 N DIHE HT2 N CA C HT1 N CA C HT3 N CA C IC HT1 N CA C 0.0000 0.00 180.00 0.00 0.0000 IC HT2 CA *N HT1 0.0000 0.00 120.00 0.00 0.0000 IC HT3 CA *N HT2 0.0000 0.00 120.00 0.00 0.0000 ! A similar objective motivates the cter patch. The lone carbonyl ! oxygen no longer exists in the zwitter ion form and is deleted in ! the patch. PRES CTER -1.00000 ! total charge of patch residue nter is -1.0 GROU ATOM C C 0.14 ATOM OT1 OC -0.57 ATOM OT2 OC -0.57 DELETE ATOM O BOND C OT1 C OT2 DIHE N CA C OT2 IMPH C CA OT2 OT1 IC N CA C OT2 0.0 0.0 180.0 0.0 0.0 IC OT2 CA *C OT1 0.0 0.0 180.0 0.0 0.0 END ! terminate the rtf reading ! Now lets print the rtf information and see what sort of connectivity ! information we get. print rtf ! information will come out in output file, not here. ! We now need to specify the parameters to be used in the potential ! energy function. These include harmonic bond terms, harmonic angle ! terms, dihedral and improper torsion terms and van der Walls constants read parameter card ! translation - read parameter file using ASCII format * - parameter file from PARAM19 - * BOND ! atom pair, bond force constant, equilibrium bond distance are given for ! each unique bond type specified in the topology. Therefore, one must ! always read a topology file first. C CH1E 405.0 1.52 C NH1 471.0 1.33 C O 580.0 1.23 C OC 580.0 1.23 CH1E CH3E 225.0 1.52 CH1E NH1 422.0 1.45 CH1E NH3 422.0 1.45 H NH1 405.0 0.98 HC NH3 405.0 1.04 ! Now all parameters for angles must be specified THETAS ! triple of atoms for angle type, force constant and equilibrium angle CH1E C NH1 20.0 117.5 CH1E C O 85.0 121.5 CH1E C OC 85.0 117.5 NH1 C O 65.0 121.0 C CH1E CH3E 70.0 106.5 C CH1E NH1 45.0 111.6 C CH1E NH3 45.0 111.6 ! We add these angle terms for the spanning angle not ! previously included. C NH1 CH1E 77.5 120.0 C NH1 H 30.0 120.0 CH1E NH1 H 35.0 120.0 ! CH3E CH1E NH1 65.0 108.5 CH3E CH1E NH3 65.0 109.5 CH1E NH3 HC 35.0 109.5 HC NH3 HC 40.0 109.5 OC C OC 85.0 122.5 PHI ! Now for the dihedrals. ! quadruple of atom types, force constant, periodicity, phase X C CH1E X 0.0 3 0.0 X C NH1 X 8.2 2 180.0 X CH1E NH1 X 0.3 3 0.0 X CH1E NH3 X 0.6 3 0.0 IMPHI ! and now impropers ! quadruple of atom types, force constant, periodicity, minimum angle C X X C 25.0 0 0.0 C X X O 100.0 0 0.0 C X X OC 100.0 0 0.0 CH1E X X CH3E 55.0 0 35.26439 ! Again, this is added to the previous parameters to make the complete for ! polypeptide. NH1 X X H 45.0 0 0.0 ! ! Lastly the non-bonded parameters, first we specify a non-bonded ! default for calculations of vdW and Electrostatic interactions ! this is not necessary as we can specify it in the specific commands ! which use the energy function. However, its a good practice to ! set up your favorite defaults. NONBONDED NBXMOD 5 ATOM CDIEL SHIFT VATOM VDISTANCE VSHIFT - CUTNB 8.0 CTOFNB 7.5 CTONNB 6.5 EPS 1.0 E14FAC 0.4 WMIN 1.5 ! ! Emin Rmin/2 ! (kcal/mol) (A) H 0.0440 -0.0498 0.8000 HC 0.0440 -0.0498 0.6000 ! charged group. Reduced vdw radius ! ! ***** Emin Rmin/2 for 1-4 interactions C 1.65 -0.1200 2.100 1.65 -0.1 1.9 ! carbonyl carbon CH1E 1.35 -0.0486 2.365 1.35 -0.1 1.9 ! \ CH3E 2.17 -0.1811 2.165 1.77 -0.1 1.9 ! / ! ! The * appended to N and O below serve as a wild card specifier which permits ! any string following the first character to be matched. Other wildcards are ! also available to match any character beyond the first. N* 1.1000 -0.2384 1.6000 ! includes N,NC2,NH1,NH2,NH3,NP,and NR ! O* 0.8400 -0.1591 1.6000 ! includes O, OH1, OM, and OS OC 2.1400 -0.6469 1.6000 END ! Now onto the science!!! We'll first build the Protein Structure File (PSF) ! This consists of a statement of the sequence, which in this case is ALA, ! followed by generation of the PSF. read sequence_of_residues card * Here we put some title again for the sequence file * Hexa-ala in a helix * 6 ! Number of residues in sequence ALA ALA ALA ALA ALA ALA ! listing of all residues in sequence. ! Note since we specified ! that the default first and last patches were to be nter and cter ! we need not specify those "residues" in the sequence. ! Here we generate the PSF and have IC tables SETUP to be used to build ! the three dimensional structure. generate hala setup ! generate a segment which we name hala ! (helical poly-alanine) and setup IC tables for ! later use. ! We can look at the PSF, it gives all the specific information regarding ! this molecular structure. But its long now so we won't! ! Now we'll write the PSF file as we did for the topology and parameters. open unit 1 write form name hala.psf write psf card unit 1 * PSF for a-helical hexa-alanine in zwitterion form * close unit 1 ! At this point we are ready to "build" a 3-D structure for the peptide from ! the information on the PSF and parameter files. We first need to fill ! all of the unspecified bond and angles in the IC tables with values ! from the parameters. The ic param command does this. ic param print ic ! Next we need to tell CHARMM's builder where to begin placing atoms. ! We'll start from the left end, with the heavy atoms. ic seed 1 N 1 CA 1 C ! specify id's of three atoms to "seed" the building ! Now build ic build coor print ! We now have a three dimensional structure for the peptide ! the coordinates can be manipulated, written out or dynamicized (?)! ! Let's write it out for posterities sake! open unit 1 write form name hala0.chr write coor card unit 1 * Initial coordinates in a-helical conformation * ! Next we will carry-out molecular dynamics for 10000 steps (10ps). ! We'll save the atomic coordinates every 100 steps to look at the ! dynamics. In general, we may want to save many more coordinate sets ! to analyze the properties of the peptide. ! We need to open the dynamics coordinate and restart units to save ! the pertinent data. The coordinate file will contain the atomic ! positions as a function of time and the restart file will contain all ! the information required to continue the calculation after quitting, i.e., ! the velocities, positions, accumulated average properties written as ! output, etc. open unit 10 write unform name hala.crd ! open a fortran unit to write a ! binary coordinate file open unit 9 write form name hala.rest ! open a formatted fortran unit to ! save the restart file ! Now do the dynamics. dynamics verlet strt nstep 10000 timestep 0.001 - ! specify dynamics control iuncrd 10 nsavc 100 iunwri 9 nprint 100 iprfrq 100 - ! output frequencies firstt 300.0 finalt 300.0 twindh 10.0 twindl -10.0 ieqfrq 50! finally, ! specify parameters for thermal control ! We can now write out the minimized coordinates with a title which ! includes the energy of this structure. open unit 1 write form name halaf.chr !hala coordinates in CHRarmm format write coor card unit 1 * a-helical hexa-alanine polypeptide zwitterion after 10 ps of MD * The final energy is ?ENER kcal/mol * stop