[CP2K:1541] ewald_type none

[Laino Teodoro]

Dear Chris,

EWALD_TYPE NONE means compute energy and forces with standard Coulomb  
potential BUT shifting the Coulomb
potential in such a way that it is zero at the RCUT_NONBONDED.
Therefore you see forces consistent but energy not (just for a  
constant). It is correct in this way.. there's no error.

Best Regards,
Teo

On 10 Nov 2008, at 20:26, Christopher Knight wrote:

>
>
> 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
>
>
>
> >