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)
          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:41:28 CST