[gmx-users] umbrella sampling Protein_ligand complex distance calculation

Discussion:

Ashwini Londhe

2018-01-12 07:14:30 UTC

Hi

Initially I followed US-gromacs tuitorial of Dr. Lemkul.it works very fine. so I following the same tutorial for ligand-protein complex.

I have pulled the ligand over 300ps using following MD_pull.mdp and also generated series of cordinate files (conf0.gro-conf300.gro) but problem is at distance calculation.

title = Umbrella pulling simulation

define = -DPOSRES

; Run parameters

integrator = md

dt = 0.002

tinit = 0

nsteps = 150000 ; 300 ps

nstcomm = 10

; Output parameters

nstxout = 5000 ; every 10 ps

nstvout = 5000

nstfout = 500

nstxtcout = 500 ; every 1 ps

nstenergy = 500

; Bond parameters

constraint_algorithm = lincs

constraints = all-bonds

continuation = yes ; continuing from NPT

; Single-range cutoff scheme

nstlist = 5

ns_type = grid

rlist = 1.4

rcoulomb = 1.4

rvdw = 1.4

; PME electrostatics parameters

coulombtype = PME

fourierspacing = 0.12

fourier_nx = 0

fourier_ny = 0

fourier_nz = 0

pme_order = 4

ewald_rtol = 1e-5

optimize_fft = yes

; Berendsen temperature coupling is on in two groups

Tcoupl = Nose-Hoover

tc_grps = Protein Non-Protein

tau_t = 0.5 0.5

ref_t = 310 310

; Pressure coupling is on

Pcoupl = Parrinello-Rahman

pcoupltype = isotropic

tau_p = 1.0

compressibility = 4.5e-5

ref_p = 1.0

refcoord_scaling = com

; Generate velocities is off

gen_vel = no

; Periodic boundary conditions are on in all directions

pbc = xyz

; Long-range dispersion correction

DispCorr = EnerPres

; Pull code

pull = yes

pull_ngroups = 2

pull_ncoords = 1

pull_group1_name = Protein

pull_group2_name = L6I

pull_coord1_type = umbrella ; harmonic biasing force

pull_coord1_geometry = distance ; simple distance increase

pull_coord1_groups = 1 2

pull_coord1_dim = N N Y

pull_coord1_rate = 0.01 ; 0.01 nm per ps = 10 nm per ns

pull_coord1_k = 1000 ; kJ mol^-1 nm^-2

pull_coord1_start = yes ; define initial COM distance > 0

　

I have used distances.pl file from gromacs tutorial and change the group name Chain_A and Chain_B with L6I and Protein. Also replaced 500 frames with 300 frames. It generated summary_distance.dat file which does not include distances values.

After running script it generates this at the end

.

.

.

Processing configuration 300...

readline () on closed filehandle IN at distances.pl line 16.

Use of uninitializd value $distance in concatenaton (.) or string at distans.pl line 30

How could it happen by just changing Chain_A to L6I and Chain_B to Protein and conformations 500 to 300 it seems script works fine for peptide system only and not for protein_ligand complex.

Also above error generated for ony proein-ligand complex and not for the peptide tutoral.

Please help me to sort out this problem.

Regards

Ashwini Londhe

Korea Institute of Science and Technology (KIST)
5, Hwarang-ro 14-gil
Seongbuk-gu
Seoul, 02792
Republic of Korea

(우: 02792) 서울특별시 성북구 화랑로 14길 5

--
Gromacs Users mailing list

* Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-***@gromacs.

rose rahmani

2018-01-12 12:09:37 UTC

Permalink

I had this problem but i couldn't solve it.
But i can tell you just a temporary solution; Don't use -sep with
trjconv,then use the output conf as an input for gmx distance. it probably
gives you sth like this;

time distances
0.0000000 1.7537287 0.0238601 0.0428935 -1.7530417
2.0000000 1.7294775 -0.0426260 -0.0664551 -1.7276745
4.0000000 1.7153677 0.0664165 0.0698769 -1.7126565
.
.
.
.
just consider that time=2 >>> is 1 ps ( depends on dt )
=4 >>> is 2ps
= 6 >>> is 3ps .........
but distances are exactly the same.
Sorrry, i know maybe Dr.justin and all professionals don't agree with these
solutions, but i can understand,when you don't get a proper answer,
sometimes you have to ridiculous solutions.That's it!

Post by Ashwini Londhe
Hi
Initially I followed US-gromacs tuitorial of Dr. Lemkul.it works very
fine. so I following the same tutorial for ligand-protein complex.
I have pulled the ligand over 300ps using following MD_pull.mdp and also
generated series of cordinate files (conf0.gro-conf300.gro) but problem is
at distance calculation.
title = Umbrella pulling simulation
define = -DPOSRES
; Run parameters
integrator = md
dt = 0.002
tinit = 0
nsteps = 150000 ; 300 ps
nstcomm = 10
; Output parameters
nstxout = 5000 ; every 10 ps
nstvout = 5000
nstfout = 500
nstxtcout = 500 ; every 1 ps
nstenergy = 500
; Bond parameters
constraint_algorithm = lincs
constraints = all-bonds
continuation = yes ; continuing from NPT
; Single-range cutoff scheme
nstlist = 5
ns_type = grid
rlist = 1.4
rcoulomb = 1.4
rvdw = 1.4
; PME electrostatics parameters
coulombtype = PME
fourierspacing = 0.12
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
pme_order = 4
ewald_rtol = 1e-5
optimize_fft = yes
; Berendsen temperature coupling is on in two groups
Tcoupl = Nose-Hoover
tc_grps = Protein Non-Protein
tau_t = 0.5 0.5
ref_t = 310 310
; Pressure coupling is on
Pcoupl = Parrinello-Rahman
pcoupltype = isotropic
tau_p = 1.0
compressibility = 4.5e-5
ref_p = 1.0
refcoord_scaling = com
; Generate velocities is off
gen_vel = no
; Periodic boundary conditions are on in all directions
pbc = xyz
; Long-range dispersion correction
DispCorr = EnerPres
; Pull code
pull = yes
pull_ngroups = 2
pull_ncoords = 1
pull_group1_name = Protein
pull_group2_name = L6I
pull_coord1_type = umbrella ; harmonic biasing force
pull_coord1_geometry = distance ; simple distance increase
pull_coord1_groups = 1 2
pull_coord1_dim = N N Y
pull_coord1_rate = 0.01 ; 0.01 nm per ps = 10 nm per ns
pull_coord1_k = 1000 ; kJ mol^-1 nm^-2
pull_coord1_start = yes ; define initial COM distance > 0
I have used distances.pl file from gromacs tutorial and change the group
name Chain_A and Chain_B with L6I and Protein. Also replaced 500 frames
with 300 frames. It generated summary_distance.dat file which does not
include distances values.
After running script it generates this at the end
.
.
.
Processing configuration 300...
readline () on closed filehandle IN at distances.pl line 16.
Use of uninitializd value $distance in concatenaton (.) or string at distans.pl line 30
How could it happen by just changing Chain_A to L6I and Chain_B to Protein
and conformations 500 to 300 it seems script works fine for peptide system
only and not for protein_ligand complex.
Also above error generated for ony proein-ligand complex and not for the peptide tutoral.
Please help me to sort out this problem.
Regards
Ashwini Londhe
Korea Institute of Science and Technology (KIST)
5, Hwarang-ro 14-gil
Seongbuk-gu
Seoul, 02792
Republic of Korea
(우: 02792) 서울특별시 성북구 화랑로 14길 5
--
Gromacs Users mailing list
* Please search the archive at http://www.gromacs.org/
Support/Mailing_Lists/GMX-Users_List before posting!
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or

Justin Lemkul

2018-01-13 00:37:34 UTC

Permalink

Post by rose rahmani
I had this problem but i couldn't solve it.
But i can tell you just a temporary solution; Don't use -sep with
trjconv,then use the output conf as an input for gmx distance. it probably
gives you sth like this;
time distances
0.0000000 1.7537287 0.0238601 0.0428935 -1.7530417
2.0000000 1.7294775 -0.0426260 -0.0664551 -1.7276745
4.0000000 1.7153677 0.0664165 0.0698769 -1.7126565
.
.
.
.
just consider that time=2 >>> is 1 ps ( depends on dt )
=4 >>> is 2ps
= 6 >>> is 3ps .........
but distances are exactly the same.
Sorrry, i know maybe Dr.justin and all professionals don't agree with these
solutions, but i can understand,when you don't get a proper answer,
sometimes you have to ridiculous solutions.That's it!

It's not ridiculous, but there is logic in either approach. My preferred
approach is to use trjconv -sep to quickly extract every frame, compile
the distances of each with consistent labels (frame numbers, not time)
and then you very easily (1) have the coordinate files already extracted
and (2) you know exactly which frames to use by frame number.

Without using my approach, one can easily analyze the entire trajectory
(perfectly fine) but then manually go back using trjconv -dump and
specify the desired frames by time. I find this to be tedious and not
easily automated.

-Justin

Post by rose rahmani

--
==================================================

Justin A. Lemkul, Ph.D.
Assistant Professor
Virginia Tech Department of Biochemistry

303 Engel Hall
340 West Campus Dr.
Blacksburg, VA 24061

***@vt.edu | (540) 231-3129
http://www.biochem.vt.edu/people/faculty/JustinLemkul.html

==================================================