- Posts: 14
Naproxen parameters
- ric.concu
- Topic Author
- Offline
- Fresh Boarder
I am writing in this forum looking for any kind of advice or tips that could help me in setting up the parameters for the parametrization of the Naproxen. I followed the tutorial on the home page however, I am really unable to find a stable configuration for this molecule. I am able to perform some simulation with the molecule alone. Although it seems working fine when I try to simulate in a box with water the simulation didn't go properly. This is my itp and how I have been defining the beads. I am using classical mdp from the martini tutorials.
Thank you in advance
[ moleculetype ]
; molname nrexcl
NAP 1
[ atoms ]
;id type resnr residu atom cgnr charge
1 Qd 1 NAP C1 1 -1
2 SC5 1 NAP C2 2 0
3 SC5 1 NAP C3 3 0
4 SC5 1 NAP C4 4 0
5 SC4 1 NAP C5 5 0
6 SC5 1 NAP C6 6 0
[bonds]
; i j funct length force.c.
1 2 1 0.25 1250
2 3 1 0.25 1250
3 4 1 0.24 1250
4 5 1 0.24 1250
5 6 1 0.25 1250
6 2 1 0.24 1250
[angles]
; i j k funct angle force.c.
1 2 3 2 127.00 25.0
2 3 4 2 179.00 25.0
3 4 5 2 60.00 25.0
4 5 6 2 120.00 25.0
5 6 2 2 120.00 25.0
6 2 3 2 120.00 25.0
1 2 6 2 149.00 25.0
Please Log in or Create an account to join the conversation.
- bart
- Away
- Admin
- Posts: 98
Cheers,
Bart
Please Log in or Create an account to join the conversation.
- ric.concu
- Topic Author
- Offline
- Fresh Boarder
- Posts: 14
thank you for your reply.
I ran a short simulation with Nap in water and the simulation seems to be okay, however, the final gro shows that the system exploded. The mdp options I am using are listed here below
integrator = md
dt = 0.0002
nsteps = 100000
nstcomm = 100
nstxout = 0
nstvout = 0
nstfout = 0
nstlog = 1000
nstenergy = 100
nstxout-compressed = 1000
compressed-x-precision = 100
cutoff-scheme = Verlet
nstlist = 10
ns_type = grid
pbc = xyz
rlist = 1.4
coulombtype = cutoff
rcoulomb_switch = 0.0
rcoulomb = 1.2
epsilon_r = 15 ; 2.5 (with polarizable water)
vdw_type = cutoff
rvdw_switch = 0.9
rvdw = 1.2
tcoupl = v-rescale
tc-grps = system
tau_t = 1.0
ref_t = 320
Pcoupl = parrinello-rahman
Pcoupltype = isotropic
tau_p = 12.0 ;parrinello-rahman is more stable with larger tau-p, DdJ, 20130422
compressibility = 3e-4
ref_p = 1.0
gen_vel = no
gen_temp = 320
gen_seed = 473529
constraints = none
constraint_algorithm = Lincs
Thank you for helping,
Riccardo
Please Log in or Create an account to join the conversation.
- riccardo
- Offline
- Platinum Boarder
Usually one first minimizes the system (with steepest descent), then run an NVT equilibration (using a Berendsen thermostat), then NPT eq (both Berendsen thermostat and barostats) and finally a run with settings similar to the ones you used (but with at least some 10 fs as timestep). Also what Bart mentioned is that , in case of instabilities, you could use a dt of 2 fs in the NVT and NPT equilibrations.
Also, you say you don't get any error message during the simulation / your simulation doesn't actually crash?
What makes you think the system has exploded in the last snapshot (your molecules' beads are very far from each other?)?
Please Log in or Create an account to join the conversation.
- ric.concu
- Topic Author
- Offline
- Fresh Boarder
- Posts: 14
I previously ran an Energy minimization with steepest descent and everytihng seems ok. Then I have been trying a NVT equilibration with Berendsen thermostat but it stops due too many LINCS warning. Due to this I tried a relaxation step (as in the Martini tutorial) but also this run miserably fails for LINCS warning. I am supposing that the system in the npt is blowing up because the the gro file looks like:
1NAP C1 11154.666 988.2201255.266 0.4971 1.6701 -4.3845
1NAP C2 21154.544 988.4531255.170 -8.4040 -1.9099 -2.0947
1NAP C3 31154.695 988.5661255.282 -5.2645 -3.8749 -3.9720
1NAP C4 41154.506 988.6561255.268 0.5445 10.6017 3.8307
1NAP C5 51154.391 988.6511255.130 5.0321 -3.1450 -0.6800
1NAP C6 61154.620 988.3591255.050 7.2961 -4.3438 9.9304
2NAP C1 71382.3181276.6201400.228 -2.6060 -0.0474 10.4751
2NAP C2 81382.0891276.7811400.218 -2.5068 -1.6939 -6.3560
Please Log in or Create an account to join the conversation.
- peterkroon
- Offline
- Gold Boarder
- Posts: 210
Please Log in or Create an account to join the conversation.
- ric.concu
- Topic Author
- Offline
- Fresh Boarder
- Posts: 14
Please Log in or Create an account to join the conversation.
- peterkroon
- Offline
- Gold Boarder
- Posts: 210
Please Log in or Create an account to join the conversation.
- ric.concu
- Topic Author
- Offline
- Fresh Boarder
- Posts: 14
Step 81, time 0.162 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 74139.070312, max 174495.281250 (between atoms 42 and 38)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
Please Log in or Create an account to join the conversation.
- peterkroon
- Offline
- Gold Boarder
- Posts: 210
Please also post your complete .top file.
Please Log in or Create an account to join the conversation.
- ric.concu
- Topic Author
- Offline
- Fresh Boarder
- Posts: 14
#include "martini_v2.2.itp"
#include "napCG_opt.itp"
#include "ions.itp"
[ system ]
NAP water
[ molecules ]
NAP 10
NA 10
W 5000
WF 500
Thank you all for your kind help
Please Log in or Create an account to join the conversation.
- peterkroon
- Offline
- Gold Boarder
- Posts: 210
Please Log in or Create an account to join the conversation.
- ric.concu
- Topic Author
- Offline
- Fresh Boarder
- Posts: 14
integrator = md
dt = 0.002
nsteps = 500000
nstxout = 5000
nstvout = 5000
nstlog = 5000
nstenergy = 300
nstxout-compressed = 250
compressed-x-grps = System
energygrps = System
nstlist = 20
ns-type = grid
rlist = 1.0
coulombtype = cut-off
rcoulomb = 1.0
rvdw = 1.0
tcoupl = Berendsen
tc-grps = System
tau-t = 0.2
ref-t = 300
Pcoupl = Berendsen
tau-p = 1.0
compressibility = 4.5e-5
ref-p = 1.0
gen-vel = yes
gen-temp = 300
gen-seed = 173529
constraints = all-bonds
[ moleculetype ]
; molname nrexcl
NAP 1
[ atoms ]
;id type resnr residu atom cgnr charge
1 Qd 1 NAP C1 1 -1
2 SC5 1 NAP C2 2 0
3 SC5 1 NAP C3 3 0
4 SC5 1 NAP C4 4 0
5 SC4 1 NAP C5 5 0
6 SC5 1 NAP C6 6 0
;[bonds]
;; i j funct length force.c.
;1 2 1 0.30 12500
;2 3 1 0.30 12500
;3 4 1 0.31 12500
;4 5 1 0.45 12500
;5 6 1 0.28 12500
;6 2 1 0.28 12500
[constraints]
; i j funct length
1 2 1 0.28
2 3 1 0.22
3 4 1 0.21
4 5 1 0.18
5 6 1 0.38
6 2 1 0.17
[angles]
; i j k funct angle force.c.
1 2 3 2 120.0 250
2 3 4 2 150.0 250
3 4 5 2 65.0 250
4 5 6 2 110.0 250
5 6 2 2 110.0 250
6 2 3 2 104.0 250
1 2 6 2 140.0 250
[dihedrals]
; i j k l funct angle force.c.
;1 2 3 4 3
;2 3 4 5 3
;3 4 5 6 3
;4 5 6 2 3
;1 2 6 5 3
Please Log in or Create an account to join the conversation.
- riccardo
- Offline
- Platinum Boarder
My guess is that you were grompping (gmx grompp) the files multiple times in the same folder trying to make your system work and started to lose track of which tpr was produced by which set of top, gro and mdp. Can you try a run (in a new folder) with the first itp (the one without constraints) and mdp's where you set 'constraints = none'?
Please Log in or Create an account to join the conversation.
- ric.concu
- Topic Author
- Offline
- Fresh Boarder
- Posts: 14
Please Log in or Create an account to join the conversation.
- riccardo
- Offline
- Platinum Boarder
Please Log in or Create an account to join the conversation.
- Pim
- Offline
- Expert Boarder
- Posts: 105
Please Log in or Create an account to join the conversation.
- ric.concu
- Topic Author
- Offline
- Fresh Boarder
- Posts: 14
Please Log in or Create an account to join the conversation.
- ric.concu
- Topic Author
- Offline
- Fresh Boarder
- Posts: 14
Please Log in or Create an account to join the conversation.
- Pim
- Offline
- Expert Boarder
- Posts: 105
Please Log in or Create an account to join the conversation.