Previous Thread
Next Thread
Print Thread
Page 2 of 4 1 2 3 4
Joined: Sep 2003
Posts: 8,605
Likes: 23
rmv Offline
Forum Member
Offline
Forum Member
Joined: Sep 2003
Posts: 8,605
Likes: 23
If your data is stored at 1.0 ps intervals (e.g. NSAVC 1000 for 1 fs timestep), changing the increment of the loop counter won't affect that, and the resulting time data would be incorrect. You'd also have to change the SKIP value when for the initial TRAJ command which sets up the trajectory reading, e.g. change it SKIP 5000 to extract points at a 5 ps interval instead of the 1 ps interval used when the file was generated.

rmv #304 04/18/05 09:54 PM
Joined: Jul 2004
Posts: 56
S
Forum Member
Offline
Forum Member
S
Joined: Jul 2004
Posts: 56
Thank Rick so much for the helpful instruction and example. Now I'm clear. I didn't specify an appropriate value for SKIP before, and therefore I didn't get the results right.

Thanks a lot.

Shirley

Joined: Nov 2005
Posts: 51
N
nnn Offline
Forum Member
Offline
Forum Member
N
Joined: Nov 2005
Posts: 51
Hi,

I would like to calculate the interaction between metal and ligand over the simulation. I have modified my script like this:
! Read in trajectories

open unit 51 read unform name "../../../6-prod/mgsys.dcd"

traj firstu 51 nunit 1 skip 5000

open write unit 21 form name inter_mg_cit.dat
write title unit 21
* time mg-ligand
*
set t 0.0 ! keep track of time, assume traj is from 0 to 750 ps

label loop
! get next coordinate set according to specifications above
traj read
! we have to update lists every time, things can move a lot in 5ps
update cutnb 12.0 ctofnb 11.0 ctonnb 8.0 -
atom fswitch cdie vdw vshift

! ?ener is total energy. other terms may also be extracted,
! see energy.doc, section Substitution, for more details
set e1 ?ener
! and protein ligand interaction
inte sele atom MMAT 154 MG end -
sele atom QMAT 155 N2 -
.or. atom QMAT 155 N1 -
.or. atom QMAT 155 N4 -
.or. atom QMAT 155 N3 -
.or. atom QMAT 155 C1 -
.or. atom QMAT 155 C2 -
.or. atom QMAT 155 O1 -
.or. atom QMAT 155 C3 -
.or. atom QMAT 155 H1 -
.or. atom QMAT 155 C12 -
.or. atom QMAT 155 O2 -
.or. atom QMAT 155 C10 -
.or. atom QMAT 155 C7 -
.or. atom QMAT 155 C6 -
.or. atom QMAT 155 H3 -
.or. atom QMAT 155 C5 -
.or. atom QMAT 155 CL1 -
.or. atom QMAT 155 C4 -
.or. atom QMAT 155 H2 -
.or. atom QMAT 155 C9 -
.or. atom QMAT 155 H4 -
.or. atom QMAT 155 C8 -
.or. atom QMAT 155 N5 -
.or. atom QMAT 155 H6 -
.or. atom QMAT 155 C11 -
.or. atom QMAT 155 H5 end

set e2 ?ener
write title unit 21
* @t @e1 @e2
*

incr t by 5.0
if t lt 1000.0 goto loop

stop

and my xxx.dat looks like this:

TIME MG-LIGAND
0.0 0 15.8467
5 15.8467 12.8259
10 12.8259 13.7943
15 13.7943 20.5529
20 20.5529 17.4542
25 17.4542 19.2578
30 19.2578 19.0469
35 19.0469 15.0701

In this case, I dont have interaction between resid 1 and resid 2, like in the example script. So which column (1 or 2), I could use?

Moreover, I got the positive value of energies (I suppose it must be negative). Is there something wrong with the script? Any suggestion is appreciated. Thanks in advance.

nnn #306 05/18/06 10:20 PM
Joined: Sep 2003
Posts: 8,605
Likes: 23
rmv Offline
Forum Member
Offline
Forum Member
Joined: Sep 2003
Posts: 8,605
Likes: 23
The example script has two INTE commands, yours has only one; that's why your values are repeated (the 1st energy value on each line is the same as the second on the previous line) and the first line has zero as the first energy. In this case you could simplify by having

write title unit 21
* @t ?ENER
*


and dropping the 'set e1' and 'set e2' commands. You only need one column of data in this case.

I don't think you should assume that INTE should be negative.


Rick Venable
computational chemist

rmv #307 05/18/06 10:52 PM
Joined: Nov 2005
Posts: 51
N
nnn Offline
Forum Member
Offline
Forum Member
N
Joined: Nov 2005
Posts: 51
Thank Prof. Rick.

Now I have one column for interaction energy.
In this case, the positive values of interaction energies means that ligand does not like to interact or form complex with metal indicating unstable system, is that right? First I assume that my ligand would strongly interact with metal ion. Any comment is very useful. Thanks.

nnn #308 05/18/06 11:44 PM
Joined: Sep 2003
Posts: 8,605
Likes: 23
rmv Offline
Forum Member
Offline
Forum Member
Joined: Sep 2003
Posts: 8,605
Likes: 23
The INTEraction energy is not necessarily highly correlated to the free energy of binding, and you may also have neglected additional interactions that may be important. Surely there are other atoms within a 5-10 A sphere around the metal besides that one ligand ...


Rick Venable
computational chemist

Joined: Jun 2004
Posts: 31
Z
zoe Offline
Forum Member
Offline
Forum Member
Z
Joined: Jun 2004
Posts: 31
Dear Prof.Nilsson,

If I want to compute the interaction energy including the images of my trajectory do i have to invoke the crystal environment to do this?

In this case, will the ?ener command give me the total interaction energy including the images' itneraction energy?

Best Regards,
Zoe

zoe #310 08/08/06 08:24 AM
Joined: Sep 2003
Posts: 4,849
Likes: 6
lennart Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 4,849
Likes: 6
The interaction energy command does not include image terms.
You could try to duplicate your system (GENE ... DUPL ..; COOR DUPL), apply COOR OPER commands to one of the copies and then do an interaction energy.


Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden
Joined: Feb 2005
Posts: 31
C
Forum Member
Offline
Forum Member
C
Joined: Feb 2005
Posts: 31
Dear Lennart,

thank you for providing the sript. it is very helpful. i am using it to calculate the interaction energies between protein and a series of ligands. i want to use these interaction energies as index of binding affinity of the ligands.

Could you please tell me how it works, namely what is the energy function of charmm INTE ENERGY command? i couldn't get a clear anwser from energy.doc where it is only indicated that terms such as bond,angle, dihedral, improper angle, van der waals, electrostatic, and hbond et al. are included in the function.

Many thanks in advance.

Cathy

Joined: Sep 2003
Posts: 4,849
Likes: 6
lennart Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 4,849
Likes: 6
The command is "INTE" (not INTE ENERGY). The energy function is basically the standard CHARMM potential energy, with some restrictions as indicated in energy.doc, eg some two-body energy terms are only computed if the first atom is present in the first selection and the second atom is in the second selection. In your case (protein + non-covalently bonded ligand) there will only be van der Waals and Coulomb contributions. Did you see any other non-zero terms printed in the output?


Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden
Page 2 of 4 1 2 3 4

Moderated by  lennart, rmv 

Link Copied to Clipboard
Powered by UBB.threads™ PHP Forum Software 7.7.5
(Release build 20201027)
Responsive Width:

PHP: 7.3.31-1~deb10u1 Page Time: 0.015s Queries: 34 (0.011s) Memory: 0.7854 MB (Peak: 0.8648 MB) Data Comp: Off Server Time: 2022-09-26 13:19:15 UTC
Valid HTML 5 and Valid CSS