- Posts: 210
Minimization of CG membrane proteins
- Robin
-
Topic Author
- Visitor
7 years 11 months ago #5576
by Robin
Minimization of CG membrane proteins was created by Robin
Dear martini users,
I'm currently working on 2 membrane receptors. I'm currently trying to minimize them in vacuum before equilibration in the membrane and got the following message:
Energy minimization has stopped, but the forces have not converged to the
requested precision Fmax < 10 (which may not be possible for your system). It
stopped because the algorithm tried to make a new step whose size was too
small, or there was no change in the energy since last step. Either way, we
regard the minimization as converged to within the available machine
precision, given your starting configuration and EM parameters.
Double precision normally gives you higher accuracy, but this is often not
needed for preparing to run molecular dynamics.
You might need to increase your constraint accuracy, or turn
off constraints altogether (set constraints = none in mdp file)
Steepest Descents converged to machine precision in 4979 steps,
but did not reach the requested Fmax < 10.
Potential Energy = -9.8869219e+03
Maximum force = 2.0205984e+02 on atom 358
Norm of force = 2.8898470e+01
I know it's not automatically an error but since the Epot doesn't converge, and an energy of 2.0205984e+02 seems pretty elevated to me, I don't have much confidence to do euilibration. Moreover, since my system (membrane + both proteins) will be simulated on the microsecond timescale, i'm afraid a bad minimized structure will lead to crashing simulations. Could you please enlighten me on the cause of this error? I'm pretty confident the atoms 358 isn't the real problem since it's not the only one to have such high forces. For more informations about my system, please find below my methods for "martinizing" the protein (dssp have been double checked) as well as my minimization.mdp file. I'm using version of gromacs 5.1.2.
#transformation to coarse grain
python martinize.py -f carma.average.pdb -o system-vaccum.top -x average_CG.pdb -ss avergage_equi.dssp -p backbone -ff elnedyn22 -cys P/CYS/121,P/CYS/198
#box definition
gmx editconf -f average_CG.pdb -d 1.0 -o average_CG.gro
#minimization
gmx grompp -f minimization.mdp -p system-vaccum.top -c average_CG.gro -o minimization-vaccum.tpr
gmx mdrun -deffnm minimization-vaccum -v
#minimization.mdp
define = -DPOSRES
integrator = steep
dt = 0.02
nsteps = 10000
nstxout = 0
nstvout = 0
nstlog = 100
nstxtcout = 100
xtc-precision = 10
rlist = 1.4
coulombtype = shift
rcoulomb = 1.2
epsilon_r = 15
vdw-type = shift
rvdw-switch = 0.9
rvdw = 1.2
tcoupl = v-rescale
tc-grps = w_ions Protein_POPC_CHOL
tau-t = 1.0 1.0
ref-t = 310 310
cutoff-scheme = group
Hope you can help,
Thank you for your attention.
I'm currently working on 2 membrane receptors. I'm currently trying to minimize them in vacuum before equilibration in the membrane and got the following message:
Energy minimization has stopped, but the forces have not converged to the
requested precision Fmax < 10 (which may not be possible for your system). It
stopped because the algorithm tried to make a new step whose size was too
small, or there was no change in the energy since last step. Either way, we
regard the minimization as converged to within the available machine
precision, given your starting configuration and EM parameters.
Double precision normally gives you higher accuracy, but this is often not
needed for preparing to run molecular dynamics.
You might need to increase your constraint accuracy, or turn
off constraints altogether (set constraints = none in mdp file)
Steepest Descents converged to machine precision in 4979 steps,
but did not reach the requested Fmax < 10.
Potential Energy = -9.8869219e+03
Maximum force = 2.0205984e+02 on atom 358
Norm of force = 2.8898470e+01
I know it's not automatically an error but since the Epot doesn't converge, and an energy of 2.0205984e+02 seems pretty elevated to me, I don't have much confidence to do euilibration. Moreover, since my system (membrane + both proteins) will be simulated on the microsecond timescale, i'm afraid a bad minimized structure will lead to crashing simulations. Could you please enlighten me on the cause of this error? I'm pretty confident the atoms 358 isn't the real problem since it's not the only one to have such high forces. For more informations about my system, please find below my methods for "martinizing" the protein (dssp have been double checked) as well as my minimization.mdp file. I'm using version of gromacs 5.1.2.
#transformation to coarse grain
python martinize.py -f carma.average.pdb -o system-vaccum.top -x average_CG.pdb -ss avergage_equi.dssp -p backbone -ff elnedyn22 -cys P/CYS/121,P/CYS/198
#box definition
gmx editconf -f average_CG.pdb -d 1.0 -o average_CG.gro
#minimization
gmx grompp -f minimization.mdp -p system-vaccum.top -c average_CG.gro -o minimization-vaccum.tpr
gmx mdrun -deffnm minimization-vaccum -v
#minimization.mdp
define = -DPOSRES
integrator = steep
dt = 0.02
nsteps = 10000
nstxout = 0
nstvout = 0
nstlog = 100
nstxtcout = 100
xtc-precision = 10
rlist = 1.4
coulombtype = shift
rcoulomb = 1.2
epsilon_r = 15
vdw-type = shift
rvdw-switch = 0.9
rvdw = 1.2
tcoupl = v-rescale
tc-grps = w_ions Protein_POPC_CHOL
tau-t = 1.0 1.0
ref-t = 310 310
cutoff-scheme = group
Hope you can help,
Thank you for your attention.
Please Log in or Create an account to join the conversation.
- peterkroon
-
- Offline
- Gold Boarder
Less
More
7 years 10 months ago #5606
by peterkroon
Replied by peterkroon on topic Minimization of CG membrane proteins
Hi Robin,
minimization != equilibration. What you did here is minimization. Try to proceed with equilibration; with a short (~10fs) timestep if nescessary; use a larger timestep for your production run.
Also, I wouldn't do the minimization/equilibration in vacuum, unless you have a very compelling reason to do so.
Other than that I don't see anything fundamentally wrong at a glance. After your production run, check the velocity distribution in your membrane and protein (separately); if that doesn't obey Boltzmann you'll need to make a third temperature coupling group (solvent, protein, membrane).
minimization != equilibration. What you did here is minimization. Try to proceed with equilibration; with a short (~10fs) timestep if nescessary; use a larger timestep for your production run.
Also, I wouldn't do the minimization/equilibration in vacuum, unless you have a very compelling reason to do so.
Other than that I don't see anything fundamentally wrong at a glance. After your production run, check the velocity distribution in your membrane and protein (separately); if that doesn't obey Boltzmann you'll need to make a third temperature coupling group (solvent, protein, membrane).
Please Log in or Create an account to join the conversation.
- Robin
-
Topic Author
- Visitor
7 years 10 months ago #5623
by Robin
Replied by Robin on topic Minimization of CG membrane proteins
Hi Peterkroon,
Thank you very much for your reply. Please don't worry, I know I was talking about minimzation in my previous post, excuse me for lacking clarity.
I already intend to run my whole simulation on 10 fs timestep since i'm afraid a 20fs timestep can crash due to clash between the 2 proteins when they start to dimerize.
I agree with you for the minimization in vacuum, I only did it since it was advised in the martini tutorial after "martinizing" the protein. Since I used backbones protein restraint, I guess the protein shouldn't change much.
Finally, thank you for the suggestion on looking at the Boltzmann distribution. I chosed to couple the proteins to the lipids since I saw that advice on the martini forums for systems similar to mine.
My protocol is: 1-minimization of the protein in vacuum
2-merging the 2 proteins at 10 nm distance in the same both_prot.gro file
3-inserting the both_prot.gro into a lipid bilayer with insane (solvatation and ionization in the same step)
4-minimizing the whole system (I got pretty much the same interruption message, with another high energy atom)
5-Equilibration of the whole system
My problem is mostly about equilibration, I tried a lot of different parameters but having a hard time to judge what the best parameters are: should I do the equilibration of the whole system in one step? Or in one NVT equilibration followed by an NPT equilibration (and should I use the .cpt from NVT for NPT, or the last coordinate .gro file is enough)?
I tried that last option, but my membrane got a lot of deformations (curvation of the membrane and diffusion of the lipid at the extremities).
I tried with ns_type = grid and pbc = xyz in production and got one protein moving out of the box. Without it my squared membrane changed into a circle during equilibration (lipids, waters and proteins seemed).
Finally I don't know if I should use the constraint-algorithm: LINCS since it is not mentioned in the tutorials files (membrane protein), but used in the cgmartini.nl/images/parameters/exampleMD...tini_v2.x_new-rf.mdp from martini .mdp examples.
About my equilibration, I think the criteria to judge if it's complete will be to check if pressure, temperature and volume have stabilized. Do you think I should check other critrias?
I'm sorry to bother with so many questions, I really want to be sure of the fiability of my equilibration before going into production.
Thank you very much for your reply. Please don't worry, I know I was talking about minimzation in my previous post, excuse me for lacking clarity.
I already intend to run my whole simulation on 10 fs timestep since i'm afraid a 20fs timestep can crash due to clash between the 2 proteins when they start to dimerize.
I agree with you for the minimization in vacuum, I only did it since it was advised in the martini tutorial after "martinizing" the protein. Since I used backbones protein restraint, I guess the protein shouldn't change much.
Finally, thank you for the suggestion on looking at the Boltzmann distribution. I chosed to couple the proteins to the lipids since I saw that advice on the martini forums for systems similar to mine.
My protocol is: 1-minimization of the protein in vacuum
2-merging the 2 proteins at 10 nm distance in the same both_prot.gro file
3-inserting the both_prot.gro into a lipid bilayer with insane (solvatation and ionization in the same step)
4-minimizing the whole system (I got pretty much the same interruption message, with another high energy atom)
5-Equilibration of the whole system
My problem is mostly about equilibration, I tried a lot of different parameters but having a hard time to judge what the best parameters are: should I do the equilibration of the whole system in one step? Or in one NVT equilibration followed by an NPT equilibration (and should I use the .cpt from NVT for NPT, or the last coordinate .gro file is enough)?
I tried that last option, but my membrane got a lot of deformations (curvation of the membrane and diffusion of the lipid at the extremities).
I tried with ns_type = grid and pbc = xyz in production and got one protein moving out of the box. Without it my squared membrane changed into a circle during equilibration (lipids, waters and proteins seemed).
Finally I don't know if I should use the constraint-algorithm: LINCS since it is not mentioned in the tutorials files (membrane protein), but used in the cgmartini.nl/images/parameters/exampleMD...tini_v2.x_new-rf.mdp from martini .mdp examples.
About my equilibration, I think the criteria to judge if it's complete will be to check if pressure, temperature and volume have stabilized. Do you think I should check other critrias?
I'm sorry to bother with so many questions, I really want to be sure of the fiability of my equilibration before going into production.
Please Log in or Create an account to join the conversation.
- peterkroon
-
- Offline
- Gold Boarder
Less
More
- Posts: 210
7 years 10 months ago #5628
by peterkroon
Replied by peterkroon on topic Minimization of CG membrane proteins
Hi Robin,
The Martini model is rather robust and resilient towards numerical instabilities, so you often 'get away' with just equilibrating in an NPT ensemble (instead of doing NVT first, and NPT after).
Personally I often do equilibration in two steps, first with Berendsen thermo- and barostat, and then with Nose-Hoover and Parrinello-Rahman (but other ensemble correct models, such as v-rescale, are also ok). I use the cpt from the first equilibration to start the second, and the cpt from the second to start the production run. This is because these contain the velocities and baro- and thermostat states in maximum precision.
Using LINCS for constraints is fine.
Word of warning: I don't work with membranes or lipids myself, so the rest could be very wrong :)
I think pbc=xyz is correct, but make sure your pressure coupling is semiisotropic (and your membrane is oriented correctly).
As for criteria to check: area per lipid and membrane thickness come to mind, but others may be applicable as well.
If in the end it turns out your system was not completely equilibrated when starting your production run, you can always discard the first n ns (-b flags in gromacs analysis tools).
HTH,
Peter
The Martini model is rather robust and resilient towards numerical instabilities, so you often 'get away' with just equilibrating in an NPT ensemble (instead of doing NVT first, and NPT after).
Personally I often do equilibration in two steps, first with Berendsen thermo- and barostat, and then with Nose-Hoover and Parrinello-Rahman (but other ensemble correct models, such as v-rescale, are also ok). I use the cpt from the first equilibration to start the second, and the cpt from the second to start the production run. This is because these contain the velocities and baro- and thermostat states in maximum precision.
Using LINCS for constraints is fine.
Word of warning: I don't work with membranes or lipids myself, so the rest could be very wrong :)
I think pbc=xyz is correct, but make sure your pressure coupling is semiisotropic (and your membrane is oriented correctly).
As for criteria to check: area per lipid and membrane thickness come to mind, but others may be applicable as well.
If in the end it turns out your system was not completely equilibrated when starting your production run, you can always discard the first n ns (-b flags in gromacs analysis tools).
HTH,
Peter
Please Log in or Create an account to join the conversation.
- Robin
-
Topic Author
- Visitor
7 years 10 months ago #5631
by Robin
Replied by Robin on topic Minimization of CG membrane proteins
Hi peter,
Thank you again for your clear advices.
I think the best option seems for my case is to use v-rescale with Berendsen during 5ns, follow by a second equilibration of the same duration with Parrinello-Rahman. Since i'm using position restraints on the protein backbone during the first run of equilibration (-DPOSRES), do you think it's logical to remove these backbone constraint during the second run of equilibration, in order to equilibrate and relax the protein?
Thank you for the tips on analysis, I'll compare my membrane parameters with published results in order to be sure of the quality of my equilibration.
Thank you again for your clear advices :)!
Robin
Thank you again for your clear advices.
I think the best option seems for my case is to use v-rescale with Berendsen during 5ns, follow by a second equilibration of the same duration with Parrinello-Rahman. Since i'm using position restraints on the protein backbone during the first run of equilibration (-DPOSRES), do you think it's logical to remove these backbone constraint during the second run of equilibration, in order to equilibrate and relax the protein?
Thank you for the tips on analysis, I'll compare my membrane parameters with published results in order to be sure of the quality of my equilibration.
Thank you again for your clear advices :)!
Robin
Please Log in or Create an account to join the conversation.
- peterkroon
-
- Offline
- Gold Boarder
Less
More
- Posts: 210
7 years 10 months ago #5633
by peterkroon
Replied by peterkroon on topic Minimization of CG membrane proteins
If you release the position restraints for the second equilibration you could get protein association there already. I would release them in the production run and discard the first few ns. The protein relaxation should be relatively fast anyway.
Please Log in or Create an account to join the conversation.
- Robin
-
Topic Author
- Visitor
7 years 10 months ago - 7 years 10 months ago #5635
by Robin
Replied by Robin on topic Minimization of CG membrane proteins
You have a point, but since I place my 2 protein at 10nm away from each other, it's very unlikely i'll see dimerization before at least 1 microsecond. Anyway, if you're right, I still can consider the second equilibration step as a production step and keep the first ns of production. Thank you for the advice :).
edit: after looking at the second 5ns of euilibration with no protein restraints, I can only see minors displacement of the protein, they still stay about 6nm apart so no problem about early dimerization. Anyway, my equilibration seems plausible, i'll look at the area per lipids and thickness of the membrane to ensure it. Thank you again for your precious help :)!
edit: after looking at the second 5ns of euilibration with no protein restraints, I can only see minors displacement of the protein, they still stay about 6nm apart so no problem about early dimerization. Anyway, my equilibration seems plausible, i'll look at the area per lipids and thickness of the membrane to ensure it. Thank you again for your precious help :)!
Last edit: 7 years 10 months ago by Robin.
Please Log in or Create an account to join the conversation.
Time to create page: 0.141 seconds