# How does NAMD calculate Electrostatic potential energy (no periodic boundary condition)?

From: Boyang Wang (pkuwangboyang_at_yahoo.com.cn)
Date: Mon Dec 26 2005 - 17:46:55 CST

Hi all, Merry Christmas!

I am now checking how NAMD calculates Electrostatic potential energy without periodic boundary condition. PME is used for perodic boundary condition.

I wrote a Fortran99 program to calculate the Electrostatic potential energy. The program is ready to run, as follows:

"
program main

character*26 o(10000,1),p(10000,1)
integer i1,i2
real*8 c1(10000,7),c2(10000,7)
real*8 es

open(file='PDB.pdb',unit=1)
close(unit=1)

c1(1,5)=-0.7
c1(2,5)=0.7

c2(1,5)=-0.7
c2(2,5)=0.7

es=0.

do 10 i1=1,2
do 20 i2=1,2

r=sqrt((c1(i1,1)-c2(i2,1))**2
& +(c1(i1,2)-c2(i2,2))**2
& +(c1(i1,3)-c2(i2,3))**2)

es=es+9.0*c1(i1,5)*
& c2(i2,5)*1.6/r

20 continue
10 continue

write(*,*) 'es= ',es,' eV'

661 format(a26,f12.3,f8.3,f8.3)
end
"

I calculate Electrostatic potential energy by E=9*10^9*q1*q2/r, q1 and q2 are two point charges, r is the distance between them. The system I calculate is two HCl molecules with coordinates specified in the PDB file. The Electrostatic potential energy I got by this calculation was about 2.4 kcal/mol, which was actually different from that given by the simulation. The 2 step simulation gave Electrostatic potential energy value about 3.0 kcal/mol. I changed the coordinates of HCl molecules, but my calculation never meets the value given by MD simulation.

The simulation was done as follows. One point worth of mention is that the 1-4 scaling factor doesn't matter so much as the system I calculate is just 2 HCl molecules.

The PDB file:

ATOM 1 Cl UNK A 0 0.000 0.000 0.000 1.00 1.00
ATOM 2 H UNK A 0 0.000 0.000 1.000 1.00 1.00
ATOM 3 Cl UNK A 0 4.000 0.000 0.000 1.00 1.00
ATOM 4 H UNK A 0 4.000 0.000 1.000 1.00 1.00

The PSF file:
PSF
1 !NTITLE
REMARKS original generated structure x-plor psf file

4 !NATOM

1 U 1 UNK Cl Cl -0.700000 35.5110
2 U 1 UNK H HA 0.700000 1.0110
3 U 1 UNK Cl Cl -0.700000 35.5110
4 U 1 UNK H HA 0.700000 1.0110

2 !NBOND: bonds
1 2 3 4

0 !NTHETA: angles

0 !NPHI: dihedrals

0 !NIMPHI: impropers

0 !NDON: donors

0 !NACC: acceptors

0 !NNB

The parameter file:

"
HA Cl 260.500 1.0000

!atom ignored epsilon Rmin/2 ignored eps,1-4 Rmin/2,1-4
Cl 0.000000 -0.227000 1.970000
HA 0.000000 -0.022000 1.320000
"

The configuration file:

"
firsttimestep 0

coordinates PDB.pdb
structure PSF.psf
paratypecharmm on
parameters st.inp

rigidbonds no
timestep 0.1
numsteps 2

##############################
temperature 270
outputname r1
binaryoutput yes

restartfreq 1000
dcdfreq 2500
xstFreq 2500
outputEnergies 1
outputPressure 1000

exclude scaled1-4
1-4scaling 1.0
margin 1.0
stepspercycle 1

useGroupPressure no
useFlexibleCell no
useConstantArea no

langevin on
langevinDamping 10
langevinHydrogen on
langevinTemp 270

switching on
switchdist 9.0
cutoff 10.0
pairlistdist 12.0
"

The output of this simulation is :

"
ETITLE: TS BOND ANGLE DIHED IMPRP ELECT VDW BOUNDARY MISC KINETIC TOTAL TEMP TOTAL2 TOTAL3 TEMPAVG

ENERGY: 0 0.0000 0.0000 0.0000 0.0000 3.0320 -0.2923 0.0000 0.0000 0.6593 3.3990 55.2935 3.3990 3.3990 55.2935

OPENING EXTENDED SYSTEM TRAJECTORY FILE
ENERGY: 1 0.0002 0.0000 0.0000 0.0000 3.0309 -0.2923 0.0000 0.0000 0.7218 3.4605 60.5336 3.4503 3.4886 60.5336

ENERGY: 2 0.0011 0.0000 0.0000 0.0000 3.0289 -0.2923 0.0000 0.0000 0.7489 3.4867 62.8131 3.4823 3.4962 62.8131
"

This problem looks simple but really puzzling to me.

Boyang

__________________________________________________
¸Ï¿ì×¢²áÑÅ»¢³¬´óÈÝÁ¿Ãâ·ÑÓÊÏä?
http://cn.mail.yahoo.com

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:41:28 CST