- Posts: 24
LINCS warning in CG energy minimisation
- acnash
-
Topic Author
- Offline
- Junior Boarder
Less
More
7 years 4 months ago #6447
by acnash
LINCS warning in CG energy minimisation was created by acnash
Hi all,
I am trying to perform a CG energy minimisation on a collagen protein. The initial structure was derived using martinize.py from our earlier studies of all-atom collagen models. For now, I am working with just the one collagen molecule (three polypeptide chains) in a fully solvated unit cell (using the PW CG water model). Unfortunately, the energy minimisation of my system is throwing up LINCS errors:
unning on 1 node with total 4 logical cores
Hardware detected:
CPU info:
Vendor: GenuineIntel
Brand: Intel(R) Core(TM) i5-4278U CPU @ 2.60GHz
SIMD instructions most likely to fit this hardware: AVX2_256
SIMD instructions selected at GROMACS compile time: AVX2_256
Reading file en_sol.tpr, VERSION 5.1.2 (double precision)
Using 1 MPI thread
Steepest Descents:
Tolerance (Fmax) = 1.00000e+03
Number of steps = 10000
Step= 0, Dmax= 1.0e-02 nm, Epot= 6.11946e+03 Fmax= 3.08223e+03, atom= 657
Step= 1, Dmax= 1.0e-02 nm, Epot= 2.51600e+03 Fmax= 2.78786e+03, atom= 2349
Step= 2, Dmax= 1.2e-02 nm, Epot= -1.70580e+03 Fmax= 7.45187e+03, atom= 3329
Step= 3, Dmax= 1.4e-02 nm, Epot= -3.46814e+03 Fmax= 2.96831e+03, atom= 4879
Step= 4, Dmax= 1.7e-02 nm, Epot= -8.43205e+03 Fmax= 2.95939e+03, atom= 4879
Step= 5, Dmax= 2.1e-02 nm, Epot= -1.37077e+04 Fmax= 7.39089e+03, atom= 4960
Step= 6, Dmax= 2.5e-02 nm, Epot= -1.58803e+04 Fmax= 4.41252e+03, atom= 4960
Step= 7, Dmax= 3.0e-02 nm, Epot= -2.02135e+04 Fmax= 6.54508e+03, atom= 3329
Step= 8, Dmax= 3.6e-02 nm, Epot= -2.32824e+04 Fmax= 8.13649e+03, atom= 5451
Step= 9, Dmax= 4.3e-02 nm, Epot= -2.60718e+04 Fmax= 5.62388e+03, atom= 5085
Step= 10, Dmax= 5.2e-02 nm, Epot= -3.06956e+04 Fmax= 3.01023e+03, atom= 1462
Step= 11, Dmax= 6.2e-02 nm, Epot= -3.86730e+04 Fmax= 1.37340e+04, atom= 4767
Step= 12, Dmax= 7.4e-02 nm, Epot= -4.11849e+04 Fmax= 7.10187e+03, atom= 1906
Step= 13, Dmax= 8.9e-02 nm, Epot= -4.56934e+04 Fmax= 8.99720e+03, atom= 5416
Step= 14, Dmax= 1.1e-01 nm, Epot= -5.01460e+04 Fmax= 6.60616e+03, atom= 1187
Step= 15, Dmax= 1.3e-01 nm, Epot= -5.46137e+04 Fmax= 6.63614e+04, atom= 1187
Step= 16, Dmax= 1.5e-01 nm, Epot= -5.59621e+04 Fmax= 3.57293e+04, atom= 5443
Step 17, time 0.425 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000000, max 0.000001 (between atoms 23307 and 23308)
bonds that rotated more than 45 degrees:
atom 1 atom 2 angle previous, current, constraint length
3321 3322 64.0 0.1100 0.1100 0.1100
Step= 17, Dmax= 1.8e-01 nm, Epot= -5.89601e+04 Fmax= 6.32789e+03, atom= 3329
Step= 19, Dmax= 1.1e-01 nm, Epot= -6.38244e+04 Fmax= 1.00835e+04, atom= 3329
Step= 20, Dmax= 1.3e-01 nm, Epot= -6.72839e+04 Fmax= 1.29215e+04, atom= 5370
Step= 21, Dmax= 1.6e-01 nm, Epot= -7.00732e+04 Fmax= 5.67008e+03, atom= 5374
Step= 22, Dmax= 1.9e-01 nm, Epot= -7.57990e+04 Fmax= 9.27826e+03, atom= 3588
Step= 23, Dmax= 2.3e-01 nm, Epot= -7.93349e+04 Fmax= 1.03771e+04, atom= 4960
Step= 24, Dmax= 2.8e-01 nm, Epot= -8.31597e+04 Fmax= 3.90494e+04, atom= 2274
Step= 25, Dmax= 3.3e-01 nm, Epot= -8.83416e+04 Fmax= 3.51779e+03, atom= 706
Step= 27, Dmax= 2.0e-01 nm, Epot= -9.51465e+04 Fmax= 3.69638e+04, atom= 3607
Step= 28, Dmax= 2.4e-01 nm, Epot= -9.95181e+04 Fmax= 8.01372e+03, atom= 3607
Step= 29, Dmax= 2.9e-01 nm, Epot= -1.05649e+05 Fmax= 4.31039e+03, atom= 5722
Step= 31, Dmax= 1.7e-01 nm, Epot= -1.09770e+05 Fmax= 8.75185e+03, atom= 3612
Step= 32, Dmax= 2.1e-01 nm, Epot= -1.11605e+05 Fmax= 1.51235e+05, atom= 3612
Step= 33, Dmax= 2.5e-01 nm, Epot= -1.15386e+05 Fmax= 2.02773e+04, atom= 5723
Step= 34, Dmax= 3.0e-01 nm, Epot= -1.18881e+05 Fmax= 4.74452e+03, atom= 6365
Step= 35, Dmax= 3.6e-01 nm, Epot= -1.22467e+05 Fmax= 4.88691e+04, atom= 6361
Step= 38, Dmax= 1.1e-01 nm, Epot= -1.24267e+05 Fmax= 5.63937e+04, atom= 6370
Step= 39, Dmax= 1.3e-01 nm, Epot= -1.26649e+05 Fmax= 4.48315e+03, atom= 6365
Step= 40, Dmax= 1.5e-01 nm, Epot= -1.31889e+05 Fmax= 5.23018e+03, atom= 3323
Step= 41, Dmax= 1.8e-01 nm, Epot= -1.35293e+05 Fmax= 8.84043e+03, atom= 3324
Step= 42, Dmax= 2.2e-01 nm, Epot= -1.39635e+05 Fmax= 9.77890e+03, atom= 1186
Step= 43, Dmax= 2.7e-01 nm, Epot= -1.42941e+05 Fmax= 3.29062e+03, atom= 1187
Step= 46, Dmax= 8.0e-02 nm, Epot= -1.46040e+05 Fmax= 3.62224e+03, atom= 3329
Step= 47, Dmax= 9.6e-02 nm, Epot= -1.48758e+05 Fmax= 5.88881e+03, atom= 1187
Step= 48, Dmax= 1.1e-01 nm, Epot= -1.50450e+05 Fmax= 2.13855e+04, atom= 1187
Step= 49, Dmax= 1.4e-01 nm, Epot= -1.51256e+05 Fmax= 5.76091e+03, atom= 5443
Step= 50, Dmax= 1.7e-01 nm, Epot= -1.53451e+05 Fmax= 2.55368e+04, atom= 1187
Step= 51, Dmax= 2.0e-01 nm, Epot= -1.54396e+05 Fmax= 4.03630e+03, atom= 1187
Step= 52, Dmax= 2.4e-01 nm, Epot= -1.57881e+05 Fmax= 6.32695e+04, atom= 1187
Step= 53, Dmax= 2.9e-01 nm, Epot= -1.57999e+05 Fmax= 6.44782e+04, atom= 5443
Step 54, time 1.35 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms nan, max 0.000001 (between atoms 15087 and 15088)
bonds that rotated more than 45 degrees:
atom 1 atom 2 angle previous, current, constraint length
Segmentation fault: 11
The first warning refers to a water molecule (the PW model is three beads - off the top of my head, I think two are virtual). The second warning, and where I see the segmentation fault, is with a HIS sidechain. Visual inspection of these two groups does not reveal anything unusual, with bond distance and bond angle fluctuating about the eq value. I've deleted that particular water molecule only to end up with some other water molecule throwing up a LINCS warning.
The MDP options for this run are:
title = Martini
;define = -DPOSRES
integrator = steep
dt = 0.01
nsteps = 10000
nstcomm = 100
emtol = 1000.0
;comm-grps =
nstxout = 1
nstvout = 0
nstfout = 0
nstlog = 1000
nstenergy = 1000
nstxout-compressed = 1000
compressed-x-precision = 100
cutoff-scheme = Verlet
nstlist = 20
ns_type = grid
pbc = xyz
verlet-buffer-tolerance = 0.005
coulombtype = PME ; reaction-field
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
rcoulomb = 1.1
epsilon_r = 2.5 ;15 ; 2.5 (with polarizable water)
epsilon_rf = 0
vdw_type = cutoff
vdw-modifier = Potential-shift-verlet
rvdw = 1.1
constraints = none
constraint_algorithm = lincs
lincs-warnangle = 45
As per instructed by the latest Martini material, I am using PME (rather than reaction-field) and a lower epsilon_r to deal with my polariazed water model.
I have decreased the integration time step (0.002, 0.001, 0.0001), increased (0.025), applied restraints to the protein first, but all attempts still result in LINCS warnings.
Any suggestions are welcome, including changes to the .mdp file; I am a CG novice.
Thanks
I am trying to perform a CG energy minimisation on a collagen protein. The initial structure was derived using martinize.py from our earlier studies of all-atom collagen models. For now, I am working with just the one collagen molecule (three polypeptide chains) in a fully solvated unit cell (using the PW CG water model). Unfortunately, the energy minimisation of my system is throwing up LINCS errors:
unning on 1 node with total 4 logical cores
Hardware detected:
CPU info:
Vendor: GenuineIntel
Brand: Intel(R) Core(TM) i5-4278U CPU @ 2.60GHz
SIMD instructions most likely to fit this hardware: AVX2_256
SIMD instructions selected at GROMACS compile time: AVX2_256
Reading file en_sol.tpr, VERSION 5.1.2 (double precision)
Using 1 MPI thread
Steepest Descents:
Tolerance (Fmax) = 1.00000e+03
Number of steps = 10000
Step= 0, Dmax= 1.0e-02 nm, Epot= 6.11946e+03 Fmax= 3.08223e+03, atom= 657
Step= 1, Dmax= 1.0e-02 nm, Epot= 2.51600e+03 Fmax= 2.78786e+03, atom= 2349
Step= 2, Dmax= 1.2e-02 nm, Epot= -1.70580e+03 Fmax= 7.45187e+03, atom= 3329
Step= 3, Dmax= 1.4e-02 nm, Epot= -3.46814e+03 Fmax= 2.96831e+03, atom= 4879
Step= 4, Dmax= 1.7e-02 nm, Epot= -8.43205e+03 Fmax= 2.95939e+03, atom= 4879
Step= 5, Dmax= 2.1e-02 nm, Epot= -1.37077e+04 Fmax= 7.39089e+03, atom= 4960
Step= 6, Dmax= 2.5e-02 nm, Epot= -1.58803e+04 Fmax= 4.41252e+03, atom= 4960
Step= 7, Dmax= 3.0e-02 nm, Epot= -2.02135e+04 Fmax= 6.54508e+03, atom= 3329
Step= 8, Dmax= 3.6e-02 nm, Epot= -2.32824e+04 Fmax= 8.13649e+03, atom= 5451
Step= 9, Dmax= 4.3e-02 nm, Epot= -2.60718e+04 Fmax= 5.62388e+03, atom= 5085
Step= 10, Dmax= 5.2e-02 nm, Epot= -3.06956e+04 Fmax= 3.01023e+03, atom= 1462
Step= 11, Dmax= 6.2e-02 nm, Epot= -3.86730e+04 Fmax= 1.37340e+04, atom= 4767
Step= 12, Dmax= 7.4e-02 nm, Epot= -4.11849e+04 Fmax= 7.10187e+03, atom= 1906
Step= 13, Dmax= 8.9e-02 nm, Epot= -4.56934e+04 Fmax= 8.99720e+03, atom= 5416
Step= 14, Dmax= 1.1e-01 nm, Epot= -5.01460e+04 Fmax= 6.60616e+03, atom= 1187
Step= 15, Dmax= 1.3e-01 nm, Epot= -5.46137e+04 Fmax= 6.63614e+04, atom= 1187
Step= 16, Dmax= 1.5e-01 nm, Epot= -5.59621e+04 Fmax= 3.57293e+04, atom= 5443
Step 17, time 0.425 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000000, max 0.000001 (between atoms 23307 and 23308)
bonds that rotated more than 45 degrees:
atom 1 atom 2 angle previous, current, constraint length
3321 3322 64.0 0.1100 0.1100 0.1100
Step= 17, Dmax= 1.8e-01 nm, Epot= -5.89601e+04 Fmax= 6.32789e+03, atom= 3329
Step= 19, Dmax= 1.1e-01 nm, Epot= -6.38244e+04 Fmax= 1.00835e+04, atom= 3329
Step= 20, Dmax= 1.3e-01 nm, Epot= -6.72839e+04 Fmax= 1.29215e+04, atom= 5370
Step= 21, Dmax= 1.6e-01 nm, Epot= -7.00732e+04 Fmax= 5.67008e+03, atom= 5374
Step= 22, Dmax= 1.9e-01 nm, Epot= -7.57990e+04 Fmax= 9.27826e+03, atom= 3588
Step= 23, Dmax= 2.3e-01 nm, Epot= -7.93349e+04 Fmax= 1.03771e+04, atom= 4960
Step= 24, Dmax= 2.8e-01 nm, Epot= -8.31597e+04 Fmax= 3.90494e+04, atom= 2274
Step= 25, Dmax= 3.3e-01 nm, Epot= -8.83416e+04 Fmax= 3.51779e+03, atom= 706
Step= 27, Dmax= 2.0e-01 nm, Epot= -9.51465e+04 Fmax= 3.69638e+04, atom= 3607
Step= 28, Dmax= 2.4e-01 nm, Epot= -9.95181e+04 Fmax= 8.01372e+03, atom= 3607
Step= 29, Dmax= 2.9e-01 nm, Epot= -1.05649e+05 Fmax= 4.31039e+03, atom= 5722
Step= 31, Dmax= 1.7e-01 nm, Epot= -1.09770e+05 Fmax= 8.75185e+03, atom= 3612
Step= 32, Dmax= 2.1e-01 nm, Epot= -1.11605e+05 Fmax= 1.51235e+05, atom= 3612
Step= 33, Dmax= 2.5e-01 nm, Epot= -1.15386e+05 Fmax= 2.02773e+04, atom= 5723
Step= 34, Dmax= 3.0e-01 nm, Epot= -1.18881e+05 Fmax= 4.74452e+03, atom= 6365
Step= 35, Dmax= 3.6e-01 nm, Epot= -1.22467e+05 Fmax= 4.88691e+04, atom= 6361
Step= 38, Dmax= 1.1e-01 nm, Epot= -1.24267e+05 Fmax= 5.63937e+04, atom= 6370
Step= 39, Dmax= 1.3e-01 nm, Epot= -1.26649e+05 Fmax= 4.48315e+03, atom= 6365
Step= 40, Dmax= 1.5e-01 nm, Epot= -1.31889e+05 Fmax= 5.23018e+03, atom= 3323
Step= 41, Dmax= 1.8e-01 nm, Epot= -1.35293e+05 Fmax= 8.84043e+03, atom= 3324
Step= 42, Dmax= 2.2e-01 nm, Epot= -1.39635e+05 Fmax= 9.77890e+03, atom= 1186
Step= 43, Dmax= 2.7e-01 nm, Epot= -1.42941e+05 Fmax= 3.29062e+03, atom= 1187
Step= 46, Dmax= 8.0e-02 nm, Epot= -1.46040e+05 Fmax= 3.62224e+03, atom= 3329
Step= 47, Dmax= 9.6e-02 nm, Epot= -1.48758e+05 Fmax= 5.88881e+03, atom= 1187
Step= 48, Dmax= 1.1e-01 nm, Epot= -1.50450e+05 Fmax= 2.13855e+04, atom= 1187
Step= 49, Dmax= 1.4e-01 nm, Epot= -1.51256e+05 Fmax= 5.76091e+03, atom= 5443
Step= 50, Dmax= 1.7e-01 nm, Epot= -1.53451e+05 Fmax= 2.55368e+04, atom= 1187
Step= 51, Dmax= 2.0e-01 nm, Epot= -1.54396e+05 Fmax= 4.03630e+03, atom= 1187
Step= 52, Dmax= 2.4e-01 nm, Epot= -1.57881e+05 Fmax= 6.32695e+04, atom= 1187
Step= 53, Dmax= 2.9e-01 nm, Epot= -1.57999e+05 Fmax= 6.44782e+04, atom= 5443
Step 54, time 1.35 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms nan, max 0.000001 (between atoms 15087 and 15088)
bonds that rotated more than 45 degrees:
atom 1 atom 2 angle previous, current, constraint length
Segmentation fault: 11
The first warning refers to a water molecule (the PW model is three beads - off the top of my head, I think two are virtual). The second warning, and where I see the segmentation fault, is with a HIS sidechain. Visual inspection of these two groups does not reveal anything unusual, with bond distance and bond angle fluctuating about the eq value. I've deleted that particular water molecule only to end up with some other water molecule throwing up a LINCS warning.
The MDP options for this run are:
title = Martini
;define = -DPOSRES
integrator = steep
dt = 0.01
nsteps = 10000
nstcomm = 100
emtol = 1000.0
;comm-grps =
nstxout = 1
nstvout = 0
nstfout = 0
nstlog = 1000
nstenergy = 1000
nstxout-compressed = 1000
compressed-x-precision = 100
cutoff-scheme = Verlet
nstlist = 20
ns_type = grid
pbc = xyz
verlet-buffer-tolerance = 0.005
coulombtype = PME ; reaction-field
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
rcoulomb = 1.1
epsilon_r = 2.5 ;15 ; 2.5 (with polarizable water)
epsilon_rf = 0
vdw_type = cutoff
vdw-modifier = Potential-shift-verlet
rvdw = 1.1
constraints = none
constraint_algorithm = lincs
lincs-warnangle = 45
As per instructed by the latest Martini material, I am using PME (rather than reaction-field) and a lower epsilon_r to deal with my polariazed water model.
I have decreased the integration time step (0.002, 0.001, 0.0001), increased (0.025), applied restraints to the protein first, but all attempts still result in LINCS warnings.
Any suggestions are welcome, including changes to the .mdp file; I am a CG novice.
Thanks
Please Log in or Create an account to join the conversation.
- peterkroon
-
- Offline
- Gold Boarder
Less
More
- Posts: 210
7 years 4 months ago #6448
by peterkroon
Replied by peterkroon on topic LINCS warning in CG energy minimisation
Hi,
changing the timestep for EM does exactly nothing ;)
Did grompp throw any warnings/notes/...?
It could be that your box is, eg, a bit high/low density causing instabilities.
Although it (particularly the last one) should not be needed, you can try adding the following to your mdp:
lincs-order=8
lincs-iter=2
Peter
changing the timestep for EM does exactly nothing ;)
Did grompp throw any warnings/notes/...?
It could be that your box is, eg, a bit high/low density causing instabilities.
Although it (particularly the last one) should not be needed, you can try adding the following to your mdp:
lincs-order=8
lincs-iter=2
Peter
Please Log in or Create an account to join the conversation.
- acnash
-
Topic Author
- Offline
- Junior Boarder
Less
More
- Posts: 24
7 years 4 months ago #6450
by acnash
Replied by acnash on topic LINCS warning in CG energy minimisation
Whoops! My session expired! So in brief!
Thanks for pointing out the timestep. I would say what on earth was I thinking, but I was/am pretty jet lagged. Excuse that mistake.
Ok, adding in the lincs-order fails to resolve the water LINCS warnings, but the reported rms now maxes out the NaN.
With regards to grompp I don't have any warnings besides a note on performance issues given PME settings if this is ran in parallel.
Finally, I have tried out a number of water densities by adjusting the vdw, from fully solvated (radius 0.15), 0.25 an 0.5. They all throw LINCS warnings on the water.
Thanks for pointing out the timestep. I would say what on earth was I thinking, but I was/am pretty jet lagged. Excuse that mistake.
Ok, adding in the lincs-order fails to resolve the water LINCS warnings, but the reported rms now maxes out the NaN.
With regards to grompp I don't have any warnings besides a note on performance issues given PME settings if this is ran in parallel.
Finally, I have tried out a number of water densities by adjusting the vdw, from fully solvated (radius 0.15), 0.25 an 0.5. They all throw LINCS warnings on the water.
Please Log in or Create an account to join the conversation.
- acnash
-
Topic Author
- Offline
- Junior Boarder
Less
More
- Posts: 24
7 years 4 months ago #6451
by acnash
Replied by acnash on topic LINCS warning in CG energy minimisation
I've also tried a minimisation in vacuum. This too results in a LINCS warning:
Step 17, time 0.425 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000000, max 0.000000 (between atoms 507 and 509)
bonds that rotated more than 45 degrees:
atom 1 atom 2 angle previous, current, constraint length
3321 3322 63.9 0.1100 0.1100 0.1100
I would like to think that it is my .mdp settings, rather than my initial model.
Step 17, time 0.425 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000000, max 0.000000 (between atoms 507 and 509)
bonds that rotated more than 45 degrees:
atom 1 atom 2 angle previous, current, constraint length
3321 3322 63.9 0.1100 0.1100 0.1100
I would like to think that it is my .mdp settings, rather than my initial model.
Please Log in or Create an account to join the conversation.
- peterkroon
-
- Offline
- Gold Boarder
Less
More
- Posts: 210
7 years 4 months ago #6452
by peterkroon
Replied by peterkroon on topic LINCS warning in CG energy minimisation
Lincs errors an sich are not a problem, especially during EM/equilibration. It *is* a problem if it causes your simulation to crash. So if an increase in lincs-order allows your EM to finish, it's resolved.
If not, how did you get your .itp files? I'm guessing you used Martinize. What was the exact command you used, and did it throw any warnings/errors/notes?
If not, how did you get your .itp files? I'm guessing you used Martinize. What was the exact command you used, and did it throw any warnings/errors/notes?
Please Log in or Create an account to join the conversation.
- acnash
-
Topic Author
- Offline
- Junior Boarder
Less
More
- Posts: 24
7 years 4 months ago #6458
by acnash
Replied by acnash on topic LINCS warning in CG energy minimisation
The .itp files were generated using the latest version of Martinize, from which no warnings or error were reported. I can't recall the exact command (I will dig it out from my notes).
I have, however, as you were right to point out, been playing with the lincs-order value. At a whim, I changed it to 18 and ran an energy minimization of the system in a vacuum. It resulted in one LINCS warning and a successful output:
writing lowest energy coordinates.
Steepest Descents converged to Fmax < 1000 in 91 steps
Potential Energy = -4.62149247646013e+04
Maximum force = 7.76276413603983e+02 on atom 6392
Norm of force = 3.92047419487426e+01
This is actually the furthest I've got with this system. I'll tinker with the lincs-order parameter further, then resolvate and try from there. Many thanks for your help so far.
I have, however, as you were right to point out, been playing with the lincs-order value. At a whim, I changed it to 18 and ran an energy minimization of the system in a vacuum. It resulted in one LINCS warning and a successful output:
writing lowest energy coordinates.
Steepest Descents converged to Fmax < 1000 in 91 steps
Potential Energy = -4.62149247646013e+04
Maximum force = 7.76276413603983e+02 on atom 6392
Norm of force = 3.92047419487426e+01
This is actually the furthest I've got with this system. I'll tinker with the lincs-order parameter further, then resolvate and try from there. Many thanks for your help so far.
Please Log in or Create an account to join the conversation.
- peterkroon
-
- Offline
- Gold Boarder
Less
More
- Posts: 210
7 years 4 months ago #6469
by peterkroon
Replied by peterkroon on topic LINCS warning in CG energy minimisation
As soon as you can finish an EM (with warnings or without, doesn't really matter), proceed to normal equilibration. If that blows up in your face, set lincs-order to 8; that should be more than enough.
Please Log in or Create an account to join the conversation.
- acnash
-
Topic Author
- Offline
- Junior Boarder
Less
More
- Posts: 24
7 years 4 months ago #6498
by acnash
Replied by acnash on topic LINCS warning in CG energy minimisation
So far:
EM of just the protein with lincs-order = 18. Threw one warning, but it did complete. Reran the EM on this output and dropped lincs-order to 8. No more lincs warnings, meet a threshold less than 1000.
Added in polarizable water, radius = 0.5. I'm not overly interested in the density as of yet, as I will be replacing much of it with further protein after I've managed to understand the EQ steps for this system. I ran an EM, it converged less than 1000 almost instantly. No warnings, lincs-order = 8.
I moved onto a very short NPT eq run just to see if my system will safely integrate forces. After 2ps, the system begins to throw LINCS errors. Interestingly, most of the warnings before Gromacs throws a segmentation fault are from water. Also, 2ps was only possible when I only request 1 process. I couldn't run this in parallel at all.
My mdp options, besides the changes required for polarizable water, are the latest .mdp file settings from the Martini website (they are below). I am going to try removing the protein altogether, and just run all my settings and steps on the retained water - for my own sanity. I would appreciate any thoughts on changes I could possibly make.
;
; STANDARD MD INPUT OPTIONS FOR MARTINI 2.x
; Updated 15 Jul 2015 by DdJ
;
; for use with GROMACS 5
; For a thorough comparison of different mdp options in combination with the Martini force field, see:
; D.H. de Jong et al., Martini straight: boosting performance using a shorter cutoff and GPUs, submitted.
title = Martini
; TIMESTEP IN MARTINI
; Most simulations are numerically stable with dt=40 fs,
; however better energy conservation is achieved using a
; 20-30 fs timestep.
; Time steps smaller than 20 fs are not required unless specifically stated in the itp file.
define = -DPOSRES
integrator = md
dt = 0.02
nsteps = 5000
nstcomm = 100
comm-grps =
nstxout = 100
nstvout = 0
nstfout = 0
nstlog = 1000
nstenergy = 1000
nstxout-compressed = 0
compressed-x-precision = 0
;compressed-x-grps =
energygrps = collagen solvent
; NEIGHBOURLIST and MARTINI
; To achieve faster simulations in combination with the Verlet-neighborlist
; scheme, Martini can be simulated with a straight cutoff. In order to
; do so, the cutoff distance is reduced 1.1 nm.
; Neighborlist length should be optimized depending on your hardware setup:
; updating ever 20 steps should be fine for classic systems, while updating
; every 30-40 steps might be better for GPU based systems.
; The Verlet neighborlist scheme will automatically choose a proper neighborlist
; length, based on a energy drift tolerance.
;
; Coulomb interactions can alternatively be treated using a reaction-field,
; giving slightly better properties.
; Please realize that electrostVatic interactions in the Martini model are
; not considered to be very accurate to begin with, especially as the
; screening in the system is set to be uniform across the system with
; a screening constant of 15. When using PME, please make sure your
; system properties are still reasonable.
;
; With the polarizable water model, the relative electrostatic screening
; (epsilon_r) should have a value of 2.5, representative of a low-dielectric
; apolar solvent. The polarizable water itself will perform the explicit screening
; in aqueous environment. In this case, the use of PME is more realistic.
cutoff-scheme = Verlet
nstlist = 20
ns_type = grid
pbc = xyz
verlet-buffer-tolerance = 0.005
coulombtype = PME ;reaction-field
rcoulomb = 1.1
epsilon_r = 2.5 ;15 ; 2.5 (with polarizable water)
epsilon_rf = 0
vdw_type = cutoff
vdw-modifier = Potential-shift-verlet
rvdw = 1.1
; MARTINI and TEMPERATURE/PRESSURE
; normal temperature and pressure coupling schemes can be used.
; It is recommended to couple individual groups in your system separately.
; Good temperature control can be achieved with the velocity rescale (V-rescale)
; thermostat using a coupling constant of the order of 1 ps. Even better
; temperature control can be achieved by reducing the temperature coupling
; constant to 0.1 ps, although with such tight coupling (approaching
; the time step) one can no longer speak of a weak-coupling scheme.
; We therefore recommend a coupling time constant of at least 0.5 ps.
; The Berendsen thermostat is less suited since it does not give
; a well described thermodynamic ensemble.
;
; Pressure can be controlled with the Parrinello-Rahman barostat,
; with a coupling constant in the range 4-8 ps and typical compressibility
; in the order of 10e-4 - 10e-5 bar-1. Note that, for equilibration purposes,
; the Berendsen barostat probably gives better results, as the Parrinello-
; Rahman is prone to oscillating behaviour. For bilayer systems the pressure
; coupling should be done semiisotropic.
tcoupl = v-rescale
tc-grps = collagen solvent
tau_t = 1.0 1.0
ref_t = 310 310
Pcoupl = berendsen ;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
refcoord_scaling = com
gen_vel = yes
gen_temp = 320
gen_seed = 473529
; MARTINI and CONSTRAINTS
; for ring systems and stiff bonds constraints are defined
; which are best handled using Lincs.
constraints = none
constraint_algorithm = lincs
lincs-warnangle = 45
lincs-order=8
lincs-iter=2
EM of just the protein with lincs-order = 18. Threw one warning, but it did complete. Reran the EM on this output and dropped lincs-order to 8. No more lincs warnings, meet a threshold less than 1000.
Added in polarizable water, radius = 0.5. I'm not overly interested in the density as of yet, as I will be replacing much of it with further protein after I've managed to understand the EQ steps for this system. I ran an EM, it converged less than 1000 almost instantly. No warnings, lincs-order = 8.
I moved onto a very short NPT eq run just to see if my system will safely integrate forces. After 2ps, the system begins to throw LINCS errors. Interestingly, most of the warnings before Gromacs throws a segmentation fault are from water. Also, 2ps was only possible when I only request 1 process. I couldn't run this in parallel at all.
My mdp options, besides the changes required for polarizable water, are the latest .mdp file settings from the Martini website (they are below). I am going to try removing the protein altogether, and just run all my settings and steps on the retained water - for my own sanity. I would appreciate any thoughts on changes I could possibly make.
;
; STANDARD MD INPUT OPTIONS FOR MARTINI 2.x
; Updated 15 Jul 2015 by DdJ
;
; for use with GROMACS 5
; For a thorough comparison of different mdp options in combination with the Martini force field, see:
; D.H. de Jong et al., Martini straight: boosting performance using a shorter cutoff and GPUs, submitted.
title = Martini
; TIMESTEP IN MARTINI
; Most simulations are numerically stable with dt=40 fs,
; however better energy conservation is achieved using a
; 20-30 fs timestep.
; Time steps smaller than 20 fs are not required unless specifically stated in the itp file.
define = -DPOSRES
integrator = md
dt = 0.02
nsteps = 5000
nstcomm = 100
comm-grps =
nstxout = 100
nstvout = 0
nstfout = 0
nstlog = 1000
nstenergy = 1000
nstxout-compressed = 0
compressed-x-precision = 0
;compressed-x-grps =
energygrps = collagen solvent
; NEIGHBOURLIST and MARTINI
; To achieve faster simulations in combination with the Verlet-neighborlist
; scheme, Martini can be simulated with a straight cutoff. In order to
; do so, the cutoff distance is reduced 1.1 nm.
; Neighborlist length should be optimized depending on your hardware setup:
; updating ever 20 steps should be fine for classic systems, while updating
; every 30-40 steps might be better for GPU based systems.
; The Verlet neighborlist scheme will automatically choose a proper neighborlist
; length, based on a energy drift tolerance.
;
; Coulomb interactions can alternatively be treated using a reaction-field,
; giving slightly better properties.
; Please realize that electrostVatic interactions in the Martini model are
; not considered to be very accurate to begin with, especially as the
; screening in the system is set to be uniform across the system with
; a screening constant of 15. When using PME, please make sure your
; system properties are still reasonable.
;
; With the polarizable water model, the relative electrostatic screening
; (epsilon_r) should have a value of 2.5, representative of a low-dielectric
; apolar solvent. The polarizable water itself will perform the explicit screening
; in aqueous environment. In this case, the use of PME is more realistic.
cutoff-scheme = Verlet
nstlist = 20
ns_type = grid
pbc = xyz
verlet-buffer-tolerance = 0.005
coulombtype = PME ;reaction-field
rcoulomb = 1.1
epsilon_r = 2.5 ;15 ; 2.5 (with polarizable water)
epsilon_rf = 0
vdw_type = cutoff
vdw-modifier = Potential-shift-verlet
rvdw = 1.1
; MARTINI and TEMPERATURE/PRESSURE
; normal temperature and pressure coupling schemes can be used.
; It is recommended to couple individual groups in your system separately.
; Good temperature control can be achieved with the velocity rescale (V-rescale)
; thermostat using a coupling constant of the order of 1 ps. Even better
; temperature control can be achieved by reducing the temperature coupling
; constant to 0.1 ps, although with such tight coupling (approaching
; the time step) one can no longer speak of a weak-coupling scheme.
; We therefore recommend a coupling time constant of at least 0.5 ps.
; The Berendsen thermostat is less suited since it does not give
; a well described thermodynamic ensemble.
;
; Pressure can be controlled with the Parrinello-Rahman barostat,
; with a coupling constant in the range 4-8 ps and typical compressibility
; in the order of 10e-4 - 10e-5 bar-1. Note that, for equilibration purposes,
; the Berendsen barostat probably gives better results, as the Parrinello-
; Rahman is prone to oscillating behaviour. For bilayer systems the pressure
; coupling should be done semiisotropic.
tcoupl = v-rescale
tc-grps = collagen solvent
tau_t = 1.0 1.0
ref_t = 310 310
Pcoupl = berendsen ;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
refcoord_scaling = com
gen_vel = yes
gen_temp = 320
gen_seed = 473529
; MARTINI and CONSTRAINTS
; for ring systems and stiff bonds constraints are defined
; which are best handled using Lincs.
constraints = none
constraint_algorithm = lincs
lincs-warnangle = 45
lincs-order=8
lincs-iter=2
Please Log in or Create an account to join the conversation.
- acnash
-
Topic Author
- Offline
- Junior Boarder
Less
More
- Posts: 24
7 years 4 months ago - 7 years 4 months ago #6500
by acnash
Replied by acnash on topic LINCS warning in CG energy minimisation
I was able to run an NPT simulation in parallel once I had removed the protein.So, again, it's down to the protein and/or protein-water mix. For the sake of completeness, I've added an image of the initial system just to give you a visual perspective of the box and protein geometry. Once I am satisfied with my knowledge of running just this system, I will be adding in further protein, to build up a unit cell of a collagen fibril, and thus under periodic images, a collagen fibre (hence why the box is much longer on one axis).
Furthermore, I noticed the PME warning from grompp:
The optimal PME mesh load for parallel simulations is below 0.5
and for highly parallel simulations between 0.25 and 0.33,
for higher performance, increase the cut-off and the PME grid spacing.
So to rule out the parallelisation of my system I ran an NPT simulation (using the .mdp option above), using PME as a serial job. This still failed. I switched to reaction-field, no grompp warning this time, but LINCS warnings followed by a crash in both a serial and parallel execution of the system.
This image is hidden for guests.
Furthermore, I noticed the PME warning from grompp:
The optimal PME mesh load for parallel simulations is below 0.5
and for highly parallel simulations between 0.25 and 0.33,
for higher performance, increase the cut-off and the PME grid spacing.
So to rule out the parallelisation of my system I ran an NPT simulation (using the .mdp option above), using PME as a serial job. This still failed. I switched to reaction-field, no grompp warning this time, but LINCS warnings followed by a crash in both a serial and parallel execution of the system.
This image is hidden for guests.
Please log in or register to see it.
Last edit: 7 years 4 months ago by acnash.
Please Log in or Create an account to join the conversation.
- peterkroon
-
- Offline
- Gold Boarder
Less
More
- Posts: 210
7 years 4 months ago #6546
by peterkroon
Replied by peterkroon on topic LINCS warning in CG energy minimisation
Ok.
So protein in vacuum is stable, and water in water is also stable.
How do you generate the initial solvated protein structure? Try using `gmx solvate` with an equilibrated water box.
How is the pressure behaving? You can try an NVT equilibration.
If the above does not work, reduce the timestep to 10 fs for equilibration purposes.
I've successfully run simulations of proteins in polarizable water, but that was always a cubic box. Maybe there's something funny there, but it'd surprise me.
PS. The PME note from grompp is just about efficiency, so it's not relevant here.
So protein in vacuum is stable, and water in water is also stable.
How do you generate the initial solvated protein structure? Try using `gmx solvate` with an equilibrated water box.
How is the pressure behaving? You can try an NVT equilibration.
If the above does not work, reduce the timestep to 10 fs for equilibration purposes.
I've successfully run simulations of proteins in polarizable water, but that was always a cubic box. Maybe there's something funny there, but it'd surprise me.
PS. The PME note from grompp is just about efficiency, so it's not relevant here.
Please Log in or Create an account to join the conversation.
- acnash
-
Topic Author
- Offline
- Junior Boarder
Less
More
- Posts: 24
7 years 4 months ago #6549
by acnash
Replied by acnash on topic LINCS warning in CG energy minimisation
Noted on the PME warning :-)
Ok, so I tried reducing the time step. Still LINCS warnings. I decided to then keep the reduced time step and switch to an NVT (maintaining the same thermostate i.e., v-scaling). The simulation lasted 1.44 ps before crashing. Monitoring the temp this is an idea of what's happening right before the lincs warnings:
......
1.080000 350.619232
1.090000 350.762909
1.100000 350.928802
1.110000 351.196930
1.120000 351.509064
1.130000 351.767883
1.140000 352.029022
1.150000 352.300507
1.160000 352.494873
1.170000 352.476929
1.180000 352.508087
1.190000 352.724670
1.200000 1928.730347
1.210000 3420.145996
1.220000 2978.011475
1.230000 2653.325928
1.240000 2692.044189
1.250000 2627.592285
1.260000 2590.594727
1.270000 2648.651611
1.280000 131403.031250
1.290000 260050.593750
1.300000 259327.906250
1.310000 256624.484375
1.320000 251619.718750
1.330000 244469.265625
1.340000 234922.468750
1.350000 223513.656250
1.360000 210495.234375
1.370000 196240.796875
1.380000 181398.640625
1.390000 166129.156250
1.400000 150897.031250
1.410000 123909.898438
1.420000 98393.601562
1.430000 87135.031250
1.440000 77253.984375
.......
So I switched to Berendsen NVT, and again with the temp set to 310:
0.660000 340.781158
0.670000 340.947998
0.680000 341.150146
0.690000 341.372650
0.700000 341.574310
0.710000 387.100677
0.720000 434.946289
0.730000 451.323883
0.740000 473.046539
0.750000 465.335815
0.760000 454.132507
0.770000 465.352386
0.780000 456.498138
0.790000 444.364197
0.800000 175053.625000
0.810000 330431.906250
0.820000 267994.531250
0.830000 208804.593750
0.840000 196390.546875
0.850000 200525.250000
0.860000 194095.921875
0.870000 188743.187500
0.880000 191966.343750
0.890000 183004.515625
0.900000 172753.984375
0.910000 180856.812500
The first warning is thrown at 0.48 ps, with a temp:
0.480000 362.135498
I understand there is no point in putting the pressure-coupling on without at least getting the temp coupling right, so I'll look into the playing with those parameters in particular. Again, any thoughts are always appreciated. Thanks
Ok, so I tried reducing the time step. Still LINCS warnings. I decided to then keep the reduced time step and switch to an NVT (maintaining the same thermostate i.e., v-scaling). The simulation lasted 1.44 ps before crashing. Monitoring the temp this is an idea of what's happening right before the lincs warnings:
......
1.080000 350.619232
1.090000 350.762909
1.100000 350.928802
1.110000 351.196930
1.120000 351.509064
1.130000 351.767883
1.140000 352.029022
1.150000 352.300507
1.160000 352.494873
1.170000 352.476929
1.180000 352.508087
1.190000 352.724670
1.200000 1928.730347
1.210000 3420.145996
1.220000 2978.011475
1.230000 2653.325928
1.240000 2692.044189
1.250000 2627.592285
1.260000 2590.594727
1.270000 2648.651611
1.280000 131403.031250
1.290000 260050.593750
1.300000 259327.906250
1.310000 256624.484375
1.320000 251619.718750
1.330000 244469.265625
1.340000 234922.468750
1.350000 223513.656250
1.360000 210495.234375
1.370000 196240.796875
1.380000 181398.640625
1.390000 166129.156250
1.400000 150897.031250
1.410000 123909.898438
1.420000 98393.601562
1.430000 87135.031250
1.440000 77253.984375
.......
So I switched to Berendsen NVT, and again with the temp set to 310:
0.660000 340.781158
0.670000 340.947998
0.680000 341.150146
0.690000 341.372650
0.700000 341.574310
0.710000 387.100677
0.720000 434.946289
0.730000 451.323883
0.740000 473.046539
0.750000 465.335815
0.760000 454.132507
0.770000 465.352386
0.780000 456.498138
0.790000 444.364197
0.800000 175053.625000
0.810000 330431.906250
0.820000 267994.531250
0.830000 208804.593750
0.840000 196390.546875
0.850000 200525.250000
0.860000 194095.921875
0.870000 188743.187500
0.880000 191966.343750
0.890000 183004.515625
0.900000 172753.984375
0.910000 180856.812500
The first warning is thrown at 0.48 ps, with a temp:
0.480000 362.135498
I understand there is no point in putting the pressure-coupling on without at least getting the temp coupling right, so I'll look into the playing with those parameters in particular. Again, any thoughts are always appreciated. Thanks
Please Log in or Create an account to join the conversation.
- acnash
-
Topic Author
- Offline
- Junior Boarder
Less
More
- Posts: 24
7 years 4 months ago #6550
by acnash
Replied by acnash on topic LINCS warning in CG energy minimisation
An interesting update, just trying to get the NVT working. I set the temp coupling to v-rescale and the tau-t to 0.1 0.1 (collagen solvent, groups).
This is actually the longest I've had the system run before the temperature rockets up. Essentially, the higher tau-t the sooner the temperate gets out of control.
;time (ps) temp (K)
5.860000 312.321808
5.880000 315.004181
5.900000 316.985901
5.920000 390.026642
5.940000 458.949707
5.960000 444.377960
5.980000 423.177948
6.000000 400.795197
6.020000 350.760956
6.040000 316.897125
6.060000 331.078522
6.080000 350.700745
6.100000 368.915405
6.120000 63132.300781
6.140000 117023.000000
This is actually the longest I've had the system run before the temperature rockets up. Essentially, the higher tau-t the sooner the temperate gets out of control.
;time (ps) temp (K)
5.860000 312.321808
5.880000 315.004181
5.900000 316.985901
5.920000 390.026642
5.940000 458.949707
5.960000 444.377960
5.980000 423.177948
6.000000 400.795197
6.020000 350.760956
6.040000 316.897125
6.060000 331.078522
6.080000 350.700745
6.100000 368.915405
6.120000 63132.300781
6.140000 117023.000000
Please Log in or Create an account to join the conversation.
- peterkroon
-
- Offline
- Gold Boarder
Less
More
- Posts: 210
7 years 4 months ago #6554
by peterkroon
Replied by peterkroon on topic LINCS warning in CG energy minimisation
That makes sense. With a shorter tau-t, the thermostat is more aggressive.
Have you tried creating the initial structure with an equilibrated water box? I personally think you're just seeing some sort of instability spinning out of control.
Have you tried creating the initial structure with an equilibrated water box? I personally think you're just seeing some sort of instability spinning out of control.
Please Log in or Create an account to join the conversation.
- acnash
-
Topic Author
- Offline
- Junior Boarder
Less
More
- Posts: 24
7 years 4 months ago #6559
by acnash
Replied by acnash on topic LINCS warning in CG energy minimisation
I solvated using the "-cs polwater.gro" file that comes with the Martini2.2 forcefield distribution. The file contains 5472 polarized water molecules, and I assume the initial water cluster is equilibrated.
This is my first time at Martini in over four years, I have never used its polarizable water model before. At the time I stuck with the single water bean (which doesn't appear to be around anymore?) Either way, an accurate representation of water interactions is what I want, given collagen-collagen association is thought to be mediated by interstitial water molecules.
I agree on the stability front. I wonder if it is worth going back to the energy minimisation stage and continue to push down the threshold. Thanks.
This is my first time at Martini in over four years, I have never used its polarizable water model before. At the time I stuck with the single water bean (which doesn't appear to be around anymore?) Either way, an accurate representation of water interactions is what I want, given collagen-collagen association is thought to be mediated by interstitial water molecules.
I agree on the stability front. I wonder if it is worth going back to the energy minimisation stage and continue to push down the threshold. Thanks.
Please Log in or Create an account to join the conversation.
- peterkroon
-
- Offline
- Gold Boarder
Less
More
- Posts: 210
7 years 4 months ago #6561
by peterkroon
Replied by peterkroon on topic LINCS warning in CG energy minimisation
The single water bead is still used extensively. PW is mainly used if you're interested in electrostatic effects, it may also be very applicable in this application ;)
More EM might help, might not. The goal for now is to make sure the simulation is going to run with a 20fs timestep. To that end you're pretty much free to change any parameter, even if it means you're going to get unphysical behaviour (since you're not going to analyse the equilibration anyway). So you can try to set the timestep even lower.
Beyond that, I'm also running out of ideas... If you can't get it to work, could you put the files somewhere I can look at them?
More EM might help, might not. The goal for now is to make sure the simulation is going to run with a 20fs timestep. To that end you're pretty much free to change any parameter, even if it means you're going to get unphysical behaviour (since you're not going to analyse the equilibration anyway). So you can try to set the timestep even lower.
Beyond that, I'm also running out of ideas... If you can't get it to work, could you put the files somewhere I can look at them?
Please Log in or Create an account to join the conversation.
- acnash
-
Topic Author
- Offline
- Junior Boarder
Less
More
- Posts: 24
7 years 4 months ago #6562
by acnash
Replied by acnash on topic LINCS warning in CG energy minimisation
Thanks for all of the suggestions so far. I am going to take a step back by removing the protein from the water and see if I can run a successful NVT with the same .mdp settings which are throwing a problem. I just want to try and rule the water out. I'll do the same for just protein on its own, although I think I've done this already (too much testing).
On a separate note, I've just had the time to download the full trajectory, of the last NVT which was throwing huge temperature values. I recorded every frame, so it's pretty big for 5 ps. There is a tiny part of the protein which is exploding! The coiled arrangement of the three helical polypeptide chains frays like a rope before some of the beads rocket off into space! Now, whether that is a result of the temperature going up or the temperature is a result of an impossible structure, is the problem I need to solve.
On a separate note, I've just had the time to download the full trajectory, of the last NVT which was throwing huge temperature values. I recorded every frame, so it's pretty big for 5 ps. There is a tiny part of the protein which is exploding! The coiled arrangement of the three helical polypeptide chains frays like a rope before some of the beads rocket off into space! Now, whether that is a result of the temperature going up or the temperature is a result of an impossible structure, is the problem I need to solve.
Please Log in or Create an account to join the conversation.
- acnash
-
Topic Author
- Offline
- Junior Boarder
Less
More
- Posts: 24
7 years 4 months ago #6563
by acnash
Replied by acnash on topic LINCS warning in CG energy minimisation
Just an idea of what the exploding protein looks like (it is a gif, load into a separate window if it doesn't load automatically). I've only displayed the backbone. From visual inspection the polypeptide chains do not appear to be starting in an impossible conformation, for example, I had a homology model of a protein once stick the backbone through the ring of a histidine residue!
This image is hidden for guests.
This image is hidden for guests.
Please log in or register to see it.
Please Log in or Create an account to join the conversation.
- peterkroon
-
- Offline
- Gold Boarder
Less
More
- Posts: 210
7 years 4 months ago #6618
by peterkroon
Replied by peterkroon on topic LINCS warning in CG energy minimisation
There's one more thing I can think of. Rotate your system such that the collagen is along the Z axis of the box, and set the pressure coupling to semiisotropic. If that also doesn't help, I'm going to need to have a look at your files, if possible.
Also note that my holiday starts on Monday, so response may be slow.
Peter
Also note that my holiday starts on Monday, so response may be slow.
Peter
Please Log in or Create an account to join the conversation.
- acnash
-
Topic Author
- Offline
- Junior Boarder
Less
More
- Posts: 24
7 years 4 months ago #6634
by acnash
Replied by acnash on topic LINCS warning in CG energy minimisation
Hi Peter,
Thanks for your help. The collagen in vacu simulation is still running using an NVT ensemble. However, this has been with an initial time step of 0.0005 ps for 24 hours (can't remember the number of steps off the top of my head), and then I've pushed it to 0.001 ps, and that too successfully ran for 24 hours. I've decided to stop there for a moment as I had a bit of a thought and I would appreciate your thoughts.
Am I right in saying that Martizine will assign 'ideal' bond angles based on the assigned secondary structure, and not that actual .pdb input layout itself? I selected -collagen which gave every backbone-backbone-backbone angle the same value. Crystal structures and atomistic simulations of collagen yields a very "unideal" coiling of the three helices. Some parts are quite flat (binding sites), whereas the telopeptide regions (end parts) are more like random coils. If the BB-BB-BB angles in my initial structure are quite far from the equilibrium angles presented in the .itp file could this contribute to immediate tension in the structure?
Either way, I've just knocked together a quick perl script to calculate the average angle from the atomistic structure and adjust the coarse grained .itp file. The only thing I won't have to accompany those new values (and they don't deviate far from the 'ideal' equilibrium value) are the bonded force constants. I once wrote a Genetic Algorithm that would return a population of force constants for bond stretch and bond angle to match the CG with those of an atomistic simulation (I measured the average and the stard deviation over a 50ps run). I might dig it up.
Thanks for your help. The collagen in vacu simulation is still running using an NVT ensemble. However, this has been with an initial time step of 0.0005 ps for 24 hours (can't remember the number of steps off the top of my head), and then I've pushed it to 0.001 ps, and that too successfully ran for 24 hours. I've decided to stop there for a moment as I had a bit of a thought and I would appreciate your thoughts.
Am I right in saying that Martizine will assign 'ideal' bond angles based on the assigned secondary structure, and not that actual .pdb input layout itself? I selected -collagen which gave every backbone-backbone-backbone angle the same value. Crystal structures and atomistic simulations of collagen yields a very "unideal" coiling of the three helices. Some parts are quite flat (binding sites), whereas the telopeptide regions (end parts) are more like random coils. If the BB-BB-BB angles in my initial structure are quite far from the equilibrium angles presented in the .itp file could this contribute to immediate tension in the structure?
Either way, I've just knocked together a quick perl script to calculate the average angle from the atomistic structure and adjust the coarse grained .itp file. The only thing I won't have to accompany those new values (and they don't deviate far from the 'ideal' equilibrium value) are the bonded force constants. I once wrote a Genetic Algorithm that would return a population of force constants for bond stretch and bond angle to match the CG with those of an atomistic simulation (I measured the average and the stard deviation over a 50ps run). I might dig it up.
Please Log in or Create an account to join the conversation.
- peterkroon
-
- Offline
- Gold Boarder
Less
More
- Posts: 210
7 years 3 months ago #6860
by peterkroon
Replied by peterkroon on topic LINCS warning in CG energy minimisation
I'm not intimately familiar with the workings of Martinize, so I'm not sure. It could explain it.
If you have the distribution of the angles from atomistic simulations you can fit a Boltzmann distribution of a harmonic potential to that to get your force constants.
If you have the distribution of the angles from atomistic simulations you can fit a Boltzmann distribution of a harmonic potential to that to get your force constants.
Please Log in or Create an account to join the conversation.
Time to create page: 0.131 seconds