ewald_type none

[Christopher Knight]

Hello everyone,

I've noticed something unexpected when calculating electrostatics
using EWALD_TYPE NONE.  

I place two atoms of opposite charge in a large box and ask what is
the potential energy due to electrostatics.  For the energy and
magnitude of the force on the particles, I expect the following in
atomic units:

q1 =  1.0e
q2 = -1.0e
r  =  5.669 Bohr

Energy = q1 * q2 / r = -0.1764
Force  = Energy / r  =  0.03111

If I ask CP2K for the potential energy and force using EWALD_TYPE
SPME, I get the following results:

Energy = -0.1730
Force  =  0.03109

in agreement with what I expect.


If I ask CP2K for the potential energy and force using EWALD_TYPE
NONE, I get the following results:

Energy = -0.02520
Force  =  0.03111

where the magnitude of the force is correct, but the energy is not.

Below is a small input file that reproduces this result.  I see the
same results with three different executables.  Two used the current
cvs source, compiled serially w/ gfortran(v4.3.2) w/ -O0 and -O3
optimization levels, lapack, and blas on a Linux-x86-64 machine.  The
other is a slightly older cvs source, compiled in parallel w/
intel(v9.1.###) and mkl on a Linux-x86-64 machine.  Is this result
reproducible on your machines?




I've managed to track this down to a call to the "potential_s" function
on line 459 in fist_nonbond_force.F.  

...
rab  = rab_com
rab2 = rab2_com
energy =  potential_s(c_spline_data,rab2,fscalar,spl_c,logger)
pot_nonbond = pot_nonbond + energy*fac
fscalar = fscalar*fac
fr(1) = fscalar*rab(1)
fr(2) = fscalar*rab(2)
fr(3) = fscalar*rab(3)
...

At this point in the code, the variables take the following values:

 rab:     0.0000000000000000        0.0000000000000000
-5.6691783986569284     
 rab2:    32.139583715798338     
 fac:    1.00000000000000000     
 energy:  -2.51989146947666810E-002
 fscalar:  -5.48832257515032094E-003


I'm having trouble reading the code for the "potential_s" function.  I
wonder if it is being passed the correct info, the output is being
handled correctly, or(probably the more likely) the input file needs
modification.  Could someone please check this?  Am I just making a
silly mistake?


Thanks in advance,

chris


==============================================

&GLOBAL
  PRINT_LEVEL LOW
  PREFERRED_FFT_LIBRARY FFTSG
  RUN_TYPE MD
&END GLOBAL

&MOTION
  &MD
    ENSEMBLE NVE
    STEPS 0
    TIMESTEP 0.5
    TEMPERATURE 300
  &END MD
&END MOTION

&FORCE_EVAL
  METHOD FIST
  
  &PRINT
    &FORCES ON
      FILENAME =force.txt
    &END FORCES
  &END
  
  &MM
    &FORCEFIELD        
      &CHARGE
        ATOM O
        CHARGE -1
      &END CHARGE

      &CHARGE
        ATOM H
        CHARGE  1
      &END CHARGE
      
      &NONBONDED
         &GENPOT
           ATOMS O O
           FUNCTION 0.0
           VARIABLES R
           RCUT 3.5
         &END GENPOT
         &GENPOT
           ATOMS O H
           FUNCTION 0.0
           VARIABLES R
           RCUT 3.5
         &END GENPOT
         &GENPOT
           ATOMS H H
           FUNCTION 0.0
           VARIABLES R
           RCUT 3.5
         &END GENPOT 
       &END NONBONDED
     &END FORCEFIELD
    
    &POISSON
      &EWALD
        EWALD_TYPE none !spme
        ALPHA .44
        GMAX 20
        O_SPLINE 6
      &END EWALD
    &END POISSON
  &END MM
  
  &SUBSYS
    &CELL
      A  30.00   0.00   0.00
      B   0.00  30.00   0.00
      C   0.00   0.00  30.00
      UNIT ANGSTROM
    &END CELL
    
    &COORD
      O  0.0    0.0    0.0
      H  0.0    0.0    3.0
    &END COORD 

  &END SUBSYS
&END FORCE_EVAL