[gmx-users] PULLING

[chris.neale at utoronto.ca]

First let me mention that I only scanned the mammoth manuscript  
snippet. Nevertheless, I'll try to address your questions.

1. You are physically able to use spc in place of tip3p. Whether or  
not that is a good idea is up to you to decide. Read the literature  
including the spc, tip3p, and other relevant forcefield papers in  
addition to any comparative studies.

2. No obvious question posed.

3. use the pull code. I don't know exactly what you want to do here...  
perhaps harmonically restrain the distance between the two mentioned  
groups to some specified value? No need to rename the water, just make  
an index group.

4. Who knows. Probably they did some preliminary simulations and  
discovered that this force constant did what they wanted. Your results  
can be entirely de-biased and so any force constant is technically  
rigorous, but some are better than others. You want it to be strong  
enough to hold your sampling near the defined displacements, but also  
weak enough that there is overlap between adjoining replicas. Probably  
time to read about the methodology here.

Chris.


Hi all.

I have to reproduce an experiment from the article "Identification of the
nucleophilic factors and the productive complex for the editing reaction by
leucyl-tRNA synthetase" (by Yohsuke Hagiwara a,b, Osamu Nureki c, Masaru
Tateno a,b,*). This is one of the steps of education I have to complete. All
actions described in this article were made with Amber9. I'm a newer in
modelling and I have only gromacs tools (v4.0.4), so I have downloaded amber
force field port-files (E.J. Sorin, S. Park). I want to know whether my mdp
files are reflexing the data I need. So this is a chapter from article:



"...1. MD simulations were performed using the AMBER 9 program. The parm99
force field was applied to the atoms of LeuRS and tRNALeu, and the TIP3P
model was used for the solvent water molecules. Electrostatic interactions
were calculated by the particle-mesh Ewald (PME) method, using 1.0 as the
dielectric constant. A cut-off of 12 ? was used to calculate the direct
space sum for PME; the electrostatic interactions beyond 12 ? were
calculated in reciprocal space by the fast Fourier transform method.

Thus, all electrostatic interactions between atoms were calculated. The
SHAKE algorithm was used to restrain the bond lengths involving hydrogen
atoms. The time step for integration was set to 1 fs. The temperature and
pres-sure were set using the Berendsen algorithm.

The initial coordinates used for the present MD simulations were from the
modelled structure previously constructed for LeuRS complexed with
valyl-tRNALeu. For the present system, we immersed the complex in a solvent
box comprising 49 587 water molecules, and used the periodic boundary
condition where the size of the unit box was 103.0 _ 138.3 _ 117.1 ?3. The
total number of atoms in the system was 165 739. In the initial (relaxation)
phase of the MD simulation, the water molecules were relaxed at 300 K for 10
ps, while the atoms of the protein and tRNALeu were positionally constrained
by a harmonic function using a force constant of 500 kcal/mol ?2. The force
constant was then reduced to 250, 125, 50, 25, 10, and 5 kcal/mol ?2 in six
MD simulations. The time consumed by each simulation was 2 ps...(then free
MD)...

2. This equilibrated system was used for further structural modeling to
investigate the hydrated structure relevant to the editing reaction, in
which 1 ns MD simulations were performed in the presence of constraints to
effectively explore the states. First, with respect to the atomic distance
between the carbonyl carbon of the substrate and the oxygen atom of the
identified nucleophilic water, a constrained MD simulation was performed to
investigate the

mechanisms of the approach of the water molecule. This simulation was
started by using a harmonic potential as the distance constraint, where the
force constant was set to 5 kcal/mol ?2. When the atomic distance was not
reduced any further, despite the presence of the harmonic potential in the
MD simulation (at about 470 ps), the force constant was increased up to 200
kcal/mol ?2; this was exploited to mimic the first phase of a nucleophilic
attack in the MD simulations, which enabled us to explore larger
conformational spaces than with the first principles MD simulations. A
similar protocol was applied in other constrained MD simulations for
efficient explorations of conformational spaces. To obtain a PMF with
respect to the rotation of the dihedral angle, C40-C30-O30-HO30 , we
employed the umbrella sampling technique, using 14 windows in the range of
50-180_. Umbrella sampling is a theoretical technique to efficiently search
for not only energetically favourable but also energetically-unfavourable
conformations in a phase space of a system, combined with a bias potential
(e.g. a harmonic potential) to overcome energy barriers. The force constants
of the  potential in those windows were set to 10.0 kcal/mol rad2.

The structures obtained using MD simulations were evaluated by exploiting
quantum mechanics/molecular mechanics (QM/MM) hybrid calculations in terms
of the stability and reactivity. For this purpose, our interface program was
used to connect QM (gamess) andMM(amber) engines..."



Questions:

1.First of all, can I use spc instead op tip3p? Some problems occured when I
tried to full the genbox with water... not enough memory (it's not true, I
made my box corresponding to the values mentioned above).



2. First part of MD (1.)

cpp                 =  /lib/cpp -traditional ;

define              =  -DFLEXIBLE;

constraints         =  hbonds

constraint_algorithm =  SHAKE

integrator          =  md

dt                  =  0.001

nsteps              =  20000

nstxout             =  2000

nstxtcout           =  2000

nstvout             =  2000

nstfout             =  2000

nstlog              =  2000

coulombtype         =  PME

rcoulomb            =  1.0

rcoulomb_switch     =  0

pme_order           =  4

optimize_fft        =  yes

fourierspacing      =  0.12

ns_type             =  grid

vdwtype             =  Cut-off

rlist               =  1.0

rvdw                =  1.4

nstlist             =  10

; Generate velocites is on at 300 K.

gen_vel             =  yes

gen_temp            =  300

gen_seed            =  173529

; Berendsen temperature coupling is on in four groups

Tcoupl              =  v-rescale

tau_t               =  0.1      0.1

tc_grps             =  protein  non-protein

ref_t               =  300      300

; Pressure coupling

Pcoupl              =  berendsen

tau_p               =  0.5

pcoupltype          =  isotropic

compressibility     =  13.5e-5

ref_p               =  3.0

; COM PULLING

; Pull type: no, umbrella, constraint or constant_force

pull                     = umbrella

pull_geometry            = position

pull_dim                 = Y Y Y

pull_group0              = protein

pull_group1              = non-protein

pull_vec1                = 0.0 0.0 0.0

pull_init1               = 0.0 0.0 0.0

pull_rate1               = 0.0

pull_k1                  = 2093



(In each of the six following simulations I will change just a number of
steps and pull value)

3. How can I perform a second part of simulation with a distance pull on
carbonyl carbon and water molecule? I have to rename this water and atom as
separate molecule? Any ideas?

  And the last one...for what reason they used this force in pulling, not
another?

Thank you, yours faithfully.