amber-developers: Problem with SHAKE when icfe != 0

From: Ilyas Yildirim <yildirim.pas.rochester.edu>
Date: Thu, 16 Mar 2006 21:50:26 -0700

Dear All,

I have done some tests using AMBER 9 when icfe != 0. I am attaching all of
the test cases with comments written in 0README file. Let me describe the
system which is being tested:

It is a 4X4 RNA helix, transforming from (GGCC)_2 -> (iGiGiCiC)_2, where G
and C are the natural bases while iG (isoguanine) and iC (isocytidine)
are the unnatural bases.

iso bases are created when the amino and carbonyl groups are transposed.
Here is a pic of cytidine and isocytidine:

         cytidine isocytidine

          C6 - C5 H41 C6 - C5 DH41
        / \ / / \ /
C1' - N1 C4-N4 ---- > C1' - N1 C4-O4
         \ / \ \ / \
          C2 - N3 H42 C2 - N3 DH42
          | |
          O2 N2
        / \ / \
     DH21 DH22 H21 H22

Here, DH21, DH22, DH41 and DH42 are dummy atoms used in order to create
the transformation in sander (or multisander). I have created the .prepi
files for the necessary residues (exactly in the same order; meaning, the
order of the atoms are the same. When there is H41 in the cytidine .prepi
file, there is DH41 in exactly the same line in the isocytidine .prepi
file).

I have created a .pdb file with water molecules for GGCC and iGiGiCiC (see
./pdb/ggcc_xleap.pdb and ./pdb/igigicic_xleap.pdb). Using the .prepi files
and extra.frcmod file located in ./residues directory, I have created
ggcc.inpcrd and ggcc.prmtop; igigicic.inpcrd and igigicic.prmtop.
ggcc.inpcrd and igigicic.inpcrd are exactly the same file.

Having created the necessary topology/cooradinate files, I have done 9
tests, and here are the results:

1. Minimizing the system with 2 CPUs worked fine, but does not work
with 4 or 8 CPUs (see ./min_cpu_2 and ./min_cpu_4).

2. I have done a simulation for 100 ps using 4 CPUs with
ntc=1,ntf=1,icfe=1,ntb=1 (also for icfe=2 case) (see
./ntc_1_ntf_1_icfe_1_ntb_1 and ./ntc_1_ntf_1_icfe_2_ntb_1). The output
files are out.0.5.p1 and out.0.5.p2 (GGCC and iGiGiCiC outputs,
respectively). And they are exactly the same, which is what is expected.
So this is good. (This is the case when SHAKE is OFF)

3. Same simulation but with ntb=2 (for icfe=1) (see
./ntc_1_ntf_1_icfe_1_ntb_2). The outputs out.0.5.p1 and out.0.5.p2 are
exactly the same, which is expected. So this is good too.

4. When SHAKE is turned on (see ./ntc_2_ntf_2_icfe_1_ntb_1), the system
blows up. I first thought that maybe if I used a .rst file rather the
inpcrd file, this might solve the problem. But it did not (see
./ntc_2_ntf_2_icfe_1_ntb_1_with_rst). I did a 2 ps simulation using the
.rst file created by ntc_1_ntf_1_icfe_1_ntb_1 simulation, and the system
still blowed up.

5. I did just a test using ntc=3 and ntf=3, and it did not work (as it is
happen if all of the bonds were constrained, but this is not possible for
parallel runs.

6. The same system (ntc=2,ntf=2,icfe=1,ntb=1) but with 1 CPU (see
ntc_2_ntf_2_icfe_1_ntb_1_CPU_1). System blows up (Did this test to see if
parallel runs can be done using 1 CPU).

------------------------------------------------------------------------
My thoughts:

The reason why it does not work when SHAKE is turned ON is because of the
DH bonds. I think SHAKE is just constraining the H-bonds, not DH-bonds.
DH-bonds are also like H-bonds, but the atom type for these dummy atoms
are DH, not H or H?.

So, for example, when one node constraints the H-bonds on -N-(H)2, the
DH-bonds on the other node, -O-(DH)2 is not constrained. I think this is
creating the whole problem. If there is a way to tell the program to
constrain the DH-bonds, too, the problem will disappear (I think). I have
renamed the atom type of DH to HD (so that the program will think
that it is an H-atom and therefore constrain that bond too), but it did
not work. I think SHAKE will only constrain the specific types of H atoms.

I would like to know how SHAKE realizes the H-bonds, if it is different
than the way I think.

Best,

-- 
  Ilyas Yildirim
  ---------------------------------------------------------------
  - Department of Chemisty       -				-
  - University of Rochester      -				-
  - Hutchison Hall, # B10        -				-
  - Rochester, NY 14627-0216     - Ph.:(585) 275 67 66 (Office)	-
  - http://www.pas.rochester.edu/~yildirim/			-
  ---------------------------------------------------------------

Received on Wed Apr 05 2006 - 23:49:38 PDT
Custom Search