CHARMM Development Project Forums User Discussion & Questions Energy terms, Constraints, Restraints, and Solvation Interaction energy of Dihederal

Hi,

I need to calculate the torsion energy of particular torsion angle (A-B-C-D) of the below molecule.

Quote:

---

H1 H3
| |
-A-B-C-D-
| |
H2 H4

---

(Note: In case the above pictorial representation does not appear properly. H1 and H2 are
connected to B and H3 and H4 are connected to C)

On using INTEraction and specifying atoms B and C in 2x(atom-selection).

Interaction sele atom segm 1 B end sele atom segm 1 C end


I get the following ouput..


SELRPN> 1 atoms have been selected out of ##
SELRPN> 1 atoms have been selected out of ##
INTER: Normal Atom Nonbondselected 0 from a total of ##
INTE ENR: Eval# ENERgy Delta-E GRMS
INTE INTERN: BONDs ANGLes UREY-b DIHEdrals IMPRopers
---------- --------- --------- --------- --------- ---------
INTE> 1 3.50397 2.41521 0.25094
INTE INTERN> 0.00165 0.00000 0.00000 3.50232 0.00000
---------- --------- --------- --------- --------- ---------


While issuing command

Interaction sele (atom segm 1 A) .or. (atom segm 1 B) end sele (atom segm 1 C) .or. (atom segm 1 D) end

I get the following output...

SELRPN> 2 atoms have been selected out of ##
SELRPN> 2 atoms have been selected out of ##
INTER: Normal Atom Nonbondselected 1 from a total of ##
INTE ENR: Eval# ENERgy Delta-E GRMS
INTE INTERN: BONDs ANGLes UREY-b DIHEdrals IMPRopers
INTE EXTERN: VDWaals ELEC HBONds ASP USER
---------- --------- --------- --------- --------- ---------
INTE> 2 8.46936 -4.96539 0.56936
INTE INTERN> 0.00165 0.00000 0.00000 3.50232 0.00000
INTE EXTERN> -0.07771 5.04311 0.00000 0.00000 0.00000
---------- --------- --------- --------- --------- ---------


Similarily any of the commands

Interaction sele (atom segm 1 H1) .or. (atom segm 1 B) end sele (atom segm 1 C) .or. (atom segm 1 H3) end

SELRPN> 2 atoms have been selected out of ##
SELRPN> 2 atoms have been selected out of ##
INTER: Normal Atom Nonbondselected 1 from a total of ##
INTE ENR: Eval# ENERgy Delta-E GRMS
INTE INTERN: BONDs ANGLes UREY-b DIHEdrals IMPRopers
INTE EXTERN: VDWaals ELEC HBONds ASP USER
---------- --------- --------- --------- --------- ---------
INTE> 3 3.96726 4.50209 0.23114
INTE INTERN> 0.00165 0.00000 0.00589 3.50232 0.00000
INTE EXTERN> -0.02358 0.48098 0.00000 0.00000 0.00000
---------- --------- --------- --------- --------- ---------


Interaction sele (atom segm 1 H1) .or. (atom segm 1 B) end sele (atom segm 1 C) .or. (atom segm 1 H4) end

SELRPN> 2 atoms have been selected out of ##
SELRPN> 2 atoms have been selected out of ##
INTER: Normal Atom Nonbondselected 1 from a total of ##
INTE ENR: Eval# ENERgy Delta-E GRMS
INTE INTERN: BONDs ANGLes UREY-b DIHEdrals IMPRopers
INTE EXTERN: VDWaals ELEC HBONds ASP USER
---------- --------- --------- --------- --------- ---------
INTE> 4 4.15722 -0.18995 0.27072
INTE INTERN> 0.00165 0.00000 0.00589 3.50232 0.00000
INTE EXTERN> -0.02593 0.67329 0.00000 0.00000 0.00000
---------- --------- --------- --------- --------- ---------




Interaction sele (atom segm 1 H2) .or. (atom segm 1 B) end sele (atom segm 1 C) .or. (atom segm 1 H3) end


SELRPN> 2 atoms have been selected out of ##
SELRPN> 2 atoms have been selected out of ##
INTER: Normal Atom Nonbondselected 1 from a total of ##
INTE ENR: Eval# ENERgy Delta-E GRMS
INTE INTERN: BONDs ANGLes UREY-b DIHEdrals IMPRopers
INTE EXTERN: VDWaals ELEC HBONds ASP USER
---------- --------- --------- --------- --------- ---------
INTE> 5 4.15722 0.00000 0.27072
INTE INTERN> 0.00165 0.00000 0.00589 3.50232 0.00000
INTE EXTERN> -0.02593 0.67329 0.00000 0.00000 0.00000
---------- --------- --------- --------- --------- ---------

Interaction sele (atom segm 1 H3) .or. (atom segm 1 B) end sele (atom segm 1 C) .or. (atom segm 1 H4) end


SELRPN> 2 atoms have been selected out of ##
SELRPN> 2 atoms have been selected out of ##
INTER: Normal Atom Nonbondselected 1 from a total of ##
INTE ENR: Eval# ENERgy Delta-E GRMS
INTE INTERN: BONDs ANGLes UREY-b DIHEdrals IMPRopers
INTE EXTERN: VDWaals ELEC HBONds ASP USER
---------- --------- --------- --------- --------- ---------
INTE> 6 3.96726 0.18995 0.23114
INTE INTERN> 0.00165 0.00000 0.00589 3.50232 0.00000
INTE EXTERN> -0.02358 0.48098 0.00000 0.00000 0.00000
---------- --------- --------- --------- --------- ---------





- In all the energy outputs, I observe that DIHEdrals energy is always 3.50232 kcal/mol.
It seems to me that the current implementation of Interaction energy calculation sums all
the dihedral angles energies of all possible dihedral angles around B-C bond.

Is there any way that I could get the dihedral energy only for A-B-C-D.

Please make me correct if I mis-interpret the Interaction command.

Best Wishes

In the CHARMM forcefield the energy associated with rotations around a particular bond in general contains all possible terms, just as you see. It is perhaps not so meaningful to look at just the A-B-C-D term; if you delete the hydrogens you may get the A-B-C-D energy.

Prof. Nilsson,

I issued following command..

delete atom sele (type H1 .or. type H2 .or. type H3 .or. type H4) end
Interaction sele atom segm 1 B end sele atom segm 1 C end

The error I get is

SELRPN> 4 atoms have been selected out of 16

Message from MAPIC: Atom numbers are changed.
DELTIC: 4 bonds deleted
DELTIC: 10 angles deleted
DELTIC: 12 dihedrals deleted

Warning from PSFSUM: The sum of charges ( -0.360000) is not an integer


***** LEVEL 0 WARNING FROM <PSFSUM> *****
***** Total charge not an integer
******************************************
BOMLEV ( 0) IS REACHED - TERMINATING. WRNLEV IS 5

It becomes clear that I have to play with BOMLEV, which I am not in favor of. Any other suggestion. Or I can play with BOMLEV.

Best Wishes,

This is a case where it is safe to reduce the bomlevel (-1). You can also set all charges to zero (before deleting the hydrogens).

I prefer to use BOMLEV (-1). Because I am just looking for particular dihedral energy and performing no real calculation. It worked..

The dihedral energy including contribution from all dihederal around B-C is 3.50232 kcal/mol and dihedral energy of only A-B-C-D is 3.50000 kcal/mol

Thanks for guidance Prof. Nilsson. Just out of curiosity, I would like to know the way one can change charges of certains atoms within CHARMM script without making changes in topology file.

Best Wishes,

scalar.doc

Prof. Nilsson,

Since you suggested putting charge to zero of hydrogens before deleting, so I issued..

scalar charge set 0 sele (type H1 .or. type H2 .or. type H3 .or. type H4) show end
delete atom sele (type H1 .or. type H2 .or. type H3 .or. type H4) end
Interaction sele atom segm 1 B end sele atom segm 1 C end

using default BOMLEV.
I got the following output...

Quote:

---

CHARMM> scalar charge set 0 sele (type H1 .or. type H2 .or. type H3 .or. type H4) show end
The following atoms are currently set:
SEGId RESId RESName .. TYPEs ..
SEGM 1 SEG H1 H2 H3 H4
SELRPN> 4 atoms have been selected out of 16

Warning from SCALAR: The sum of charges ( -0.360000) is not an integer


CHARMM> delete atom sele (type H1 .or. type H2 .or. type H3 .or. type H4) end
SELRPN> 4 atoms have been selected out of 16

Message from MAPIC: Atom numbers are changed.
DELTIC: 4 bonds deleted
DELTIC: 10 angles deleted
DELTIC: 12 dihedrals deleted

Warning from PSFSUM: The sum of charges ( -0.360000) is not an integer


***** LEVEL 0 WARNING FROM <PSFSUM> *****
***** Total charge not an integer
******************************************
BOMLEV ( 0) IS REACHED - TERMINATING. WRNLEV IS 5

---

It gives me an impression that even assigning charge zero to hydrogen atoms makes no difference and one has to go for BOMLEV (-1)

Best Wishes,

?
The "problem" you have is that CHARMM warns you that the total charge is non-integer, AFTER you delete the hydrogens, which is in this case not surprising and also not a real problem.
Setting the hydrogen charges to zero is not going to change this in any way since the hydrogens are already gone when CHARMM detects a possible problem with non-integer total charge.
You need to zero ALL charges for this approach to work.

issuing..

scalar charge set 0 sele all show end


solves the problem..
Thanks

best wishes,