Dear gromacs users,
I have been got stuck for a long time with a membrane system solvated in polarizable martini water model. I understand this is gromacs forum, but if anyone is aware of this issue and help me that would be really great. Any help is much appreciated.
The polarizable martini water model (martini_v2.2refP.itp) requires a different itp file compared to the normal cg martini bead model.
I am using polarizable water model. I am adding the water molecules using the following command:
gmx solvate -cp vac_em.gro -cs polarize-water.gro -p system.top -radius 0.21 -o sol.gro
Then I am adding counterions using genion command followed by energy minimization using steepest descents. I have commented the constraints section for polarizable water in the itp file and using stiff bonds with high force constants (default value provided on martini website is 50000 kJ mol-1 nm-2); However, I am not able to minimize the system, the Fmax value is in the order of 10^4 to 10^7. The energy minimization gro file I checked in vmd, that looks like the water molecules are forming some kind of clusters, almost 5-10 water molecules are almost forming too close to each other and looks very unphysical. I tried with decreasing the force constant value to 50-5000, now I get good structure of water but Fmax is super high and when I am trying equilibration the simulation is crashing. Just to check I tried the energy minimization on the polarizable waterbox provided on martini website, there also I am facing the same issue. A little less problem if I increase the force constant of the stiff bonds to 50000000 or something. I have been spending a lot of time on this without any luck. I sincerely ask for help, any help will be much appreciated, thanks in advance.
I am using the following itp file for polarizable water during the energy minimization: (martini_v2.2refP.itp)
;;;;;; POLARIZABLE WATER
[ moleculetype ]
; molname nrexcl
PW 1
[ atoms ]
;id type resnr residu atom cgnr charge
1 POL 1 PW W 1 0
2 D 1 PW WP 1 0.457
3 D 1 PW WM 1 -0.457
;[constraints]
;; i j funct length
; 1 2 1 0.14
; 1 3 1 0.14
; for minimization purposes constraints might be replaced by stiff bonds:
[bonds]
; i j funct length force const.
1 2 1 0.14 50000
1 3 1 0.14 50000
;#endif
[angles]
; i j k funct angle fc
2 1 3 2 0.0 4.20
[exclusions]
1 2 3
2 3
;;;;;; ANTIFREEZE (prevents freezing of water)
[ moleculetype ]
; molname nrexcl
WF 1
[ atoms ]
;id type resnr residu atom cgnr charge
1 BP4 1 WF WF 1 0
The energy minimization mdp file I using:
integrator = steep
;emtol = 100
nsteps = 500000
nstxout = 0
nstfout = 0
nstlog = 100
cutoff-scheme = Verlet
nstlist = 10
ns_type = grid
pbc = xyz
verlet-buffer-tolerance = 0.005
coulombtype = reaction-field
rcoulomb = 1.1
epsilon_r = 2.5 ; 2.5 (with polarizable water)
epsilon_rf = 0
vdw_type = cutoff
vdw-modifier = Potential-shift-verlet
rvdw = 1.1