1. General topics
  3. Other
  5. My Polarizable Water (PW) beads are bending after EM

normal My Polarizable Water (PW) beads are bending after EM

  • jimmy_chang
  • jimmy_chang's Avatar Topic Author
  • Offline
  • Fresh Boarder
More
6 years 3 months ago #7471 by jimmy_chang
My Polarizable Water (PW) beads are bending after EM was created by jimmy_chang
Hi there.

I seem to have a little bit of a problem during my EM process, and I really need some help.
Currently, I am trying to observe the self-assembly tendency of hydrophobin proteins. Here's my MARINI simulation process:

1. To start off I first martinized a monomer protein using the following command:
python martinize.py -f Protein.pdb -ss Protein.dssp -ff martini22p -p backbone -o topol.top -x Protein_CG.pdb
(It should be noted that my original all-atom monomer protein had 70 residues consisting 1010 atoms. After martinizing, it was reduced to 171 atoms). martinize.py was version 2.6

2. Then, I added my .itp files in my topology, which were: "martini_v2.2refP.itp," "martini_v2.2p_aminoacids," and "martini_v2.0_ions.itp" (The .itp files were from www.cgmartini.nl/index.php/force-field-p...particle-definitions
www.cgmartini.nl/index.php/force-field-parameters/amino-acids
and www.cgmartini.nl/index.php/force-field-parameters/ions )

3. I created my .gro file:
gmx editconf -f Protein_CG.pdb -c -d 1.5 -bt cubic -o Protein_CG.gro
The resulting box dimensions here were: [6.13912 6.13912 6.13912].

4. I ran a vacuum EM on this monomer:
gmx grompp -f VacuumEM.mdp -c Protein_CG.gro -p topol.top -o VacuumEM.tpr -maxwarn 1
The flag 'maxwarn 1' was for the atom name mismatch on SCD-SCN and SCD-SCP. My VacuumEM.mdp is:

Warning: Spoiler! [ Click to expand ]


You can see that my lincs-order and lincs-iter is pretty high, but I did not want to have any LINCS error during this process so I raised the value a bit. Anyways, from here I got my VacuumEM.gro file

5. After that I replicated my energy minimized structure using genconf:
gmx genconf -f VacuumEM.gro -nbox 4 4 3 -o 48-Mer.gro
This created a new box with dimensions [24.55648 24.55648 18.41736].

6. Then I solvated my system with polarizable water. First I used the following command to solvate with non-polarizable water:
gmx solvate -cp 48-Mer.gro -cs water.gro -p topol.top -o solv.gro -radius 0.21
and then ran python triple-w.py solv.gro to convert non-polarizable water to polar waters. This gave me the result gro file solv_PW.gro.
(I got water.gro and triple-w.py from www.cgmartini.nl/index.php/example-appli...ons2/solvent-systems and www.cgmartini.nl/index.php/tools2/resolution-transformation )

7. Since my system had no net charge, I skipped my ionization process and began running EM. In prior to creating a .tpr file, I read a post saying that the minimization process works better to use stiff bonds on PW molecules rather than constraints. So I have I have commented the [constraints] section and uncommented the [bonds] section:

Warning: Spoiler! [ Click to expand ]


After that, I ran an EM using 12 core workstation:
gmx grompp -f EM.mdp -c solv_PW.gro -p topol.top -o EM.tpr
mpirun -np 12 -host node01 mdrun_mpi -v -deffnm EM

Here is the EM.mdp file:

Warning: Spoiler! [ Click to expand ]


It's almost the same as VacuumEM.mdp, but I changed the 'periodic-molecules' as 'yes' along with emstep of 0.002

I anticipated my EM simulation to last a bit longer, but it finished very quickly with the following results:

Warning: Spoiler! [ Click to expand ]


8. I changed my PW stiff bonds to constraints again:
Warning: Spoiler! [ Click to expand ]


and tried to run my EQ with 10fs first.
gmx grompp -f EQ.mdp -c EM.gro -p topol.top -o EQ10.tpr
mpirun -np 12 -host node02 mdrun_mpi -v -deffnm EQ10

This is my EQ.mdp:

Warning: Spoiler! [ Click to expand ]


However, my system crashes with the following error:

Warning: Spoiler! [ Click to expand ]


Now, when I checked these LINCS erred molecules from EM.gro, I noticed that the water molecules were bent during EM

Does anyone know why this is happening? I looked up various things but I just can't seem to find out the solution.

Oh FYI, here are my topology and Protein_A.itp files:

topol.top
Warning: Spoiler! [ Click to expand ]


Protein_A.itp
Warning: Spoiler! [ Click to expand ]


Thank you very much in advance!

P.S I really wanted to upload the water molecules but I can't seem to find a way... if any knows how to upload an image please let me know!

Please Log in or Create an account to join the conversation.

More
6 years 3 months ago #7479 by Pim
Replied by Pim on topic My Polarizable Water (PW) beads are bending after EM
Hi Jimmy, thanks for your very clear and well-documented question.

You can see that near the end of your energy minimization the energy goes up steeply (to 10e12), which indeed indicates something is wrong with waters going crazy. Your NVT run therefore immediately explodes. Sadly, even with stiff bonds instead of constraints this sometimes happens. With your systems size, aiming for a max force of 10e3 or 10e4 is normal. The simplest solution that will probably work is simply taking the structure from step 112 or so and continuing with that (if you write every frame with nstxout-compressed=1 you can get it out using trjconv -dump or something like that). However, do also check if these waters are perhaps at the boundaries of your genconf-generated boxes and if the density is homogeneous.

If the above doesn't work, try solvating your system with a different random seed. In some of my self-assembly simulations, I find that in about 10% of the cases I have similar problems and just trying again works. Also, make sure you are using a gromacs version newer than 5.1, gmx solvate had bugs in it. Still even with the latest gmx sometimes I find gmx solvate doesn't always work well for this kind of setups and for some mystery reason you get a box with correct density of water at the edges, but more density in a cube in the middle of the box. You can easily spot this using VMD though. However, with your way of doing it (solvating 1 protein and then using genconf to replicate it), this problem should be minimal. I usually solvate all proteins at the same time (gmx insert-molecules the protein -nmol 48 so they are randomly rotated and then gmx solvate)

Some other tips:
- You never need to do 10,000,000 EM steps. Usually 10,000 is already more than enough.
- You don't have periodic molecules, so don't set that to yes. periodic molecules are molecules that covalently continue on both sides of the box, like infinite polymers or surfaces.
- You probably want to generate some velocities in your NVT start-up (gen_vel yes, otherwise you start at 0 Kelvin)

Please Log in or Create an account to join the conversation.

More
6 years 3 months ago #7482 by peterkroon
Replied by peterkroon on topic My Polarizable Water (PW) beads are bending after EM
In addition, did you visually check your protein after the vacuum minimization? If that ran for 10m steps it could produce an unstable conformation.
If what Pim said doesn't help put position restraints on your protein, and equilibrate your solvent (NVT or NPT doesn't matter) with a lower than normal timestep (e.g. 10 fs); or even with the stiff bonds instead of constraints.

Please Log in or Create an account to join the conversation.

  1. General topics
  3. Other
  5. My Polarizable Water (PW) beads are bending after EM
Time to create page: 0.115 seconds