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)
read(1,661) o(1,1),c1(1,1),c1(1,2),c1(1,3)
read(1,661) o(1,2),c1(2,1),c1(2,2),c1(2,3)
read(1,661) p(1,1),c2(1,1),c2(1,2),c2(1,3)
read(1,661) p(1,2),c2(2,1),c2(2,2),c2(2,3)
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
"
Added to the Charmm forcefield "par_all22_prot.inp" file, which has "HA" already.
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.
Thanks for all your time!
Boyang
__________________________________________________
¸Ï¿ì×¢²áÑÅ»¢³¬´óÈÝÁ¿Ãâ·ÑÓÊÏä?
http://cn.mail.yahoo.com
This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:40:14 CST