mirror of
https://github.com/cp2k/cp2k.git
synced 2026-07-21 14:35:15 -04:00
101 lines
3.3 KiB
Fortran
101 lines
3.3 KiB
Fortran
!-----------------------------------------------------------------------------!
|
|
! CP2K: A general program to perform molecular dynamics simulations !
|
|
! Copyright (C) 2005 CP2K developers group !
|
|
!-----------------------------------------------------------------------------!
|
|
|
|
!!****h* cp2k/f_diff_low [1.0] *
|
|
!!
|
|
!! NAME
|
|
!! f_diff_low
|
|
!!
|
|
!! FUNCTION
|
|
!! This an example of a program that uses the f77 interface.
|
|
!! This program reads the file "input.inp" that should a valid
|
|
!! cp2k input file for something that calculates the force
|
|
!! (MD for example).
|
|
!! Then this program calculates the forces with finite differences and
|
|
!! compares them to the analytic ones. The error and the percentual error
|
|
!! with respect to the mean quadratic force are printed.
|
|
!!
|
|
!! NOTES
|
|
!! You should be able to link it against the cp2k library
|
|
!! Something like
|
|
!! cd ../../makefiles
|
|
!! make libs VERSION=sopt
|
|
!! cd ../tools/f77_interface
|
|
!! f95 f_diff_low.f90 -L${HOME}/cp2k/lib/Linux-i686-nag/sdbg \
|
|
!! -lcp2k_lib -lcp2k_base_lib -L${HOME}/lib -llapack -lg2c
|
|
!! should work
|
|
!!
|
|
!! AUTHOR
|
|
!! Fawzi Mohamed
|
|
!!
|
|
!! MODIFICATION HISTORY
|
|
!! 05.2005 created [fawzi]
|
|
!!
|
|
!!*** ***********************************************************************
|
|
PROGRAM f_diff_low
|
|
IMPLICIT NONE
|
|
INTEGER :: ierr, f_env_id, natom, i
|
|
INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND ( 14, 200 )
|
|
REAL(dp), DIMENSION(:), POINTER :: pos,f
|
|
REAL(dp) :: e0,ep,em,p0,f_m,err_max,err
|
|
REAL(dp), PARAMETER :: small=5.e-3
|
|
|
|
CALL cp_init_cp2k(1,ierr)
|
|
IF (ierr/=0) STOP "init_cp2k"
|
|
|
|
CALL cp_create_fenv(f_env_id,"input.inp","out.out",ierr)
|
|
IF (ierr/=0) STOP "create_force_env"
|
|
CALL cp_get_natom(f_env_id,natom,ierr)
|
|
IF (ierr/=0) STOP "get_natom"
|
|
|
|
ALLOCATE(pos(natom*3),f(natom*3),stat=ierr)
|
|
IF (ierr/=0) STOP "alloc"
|
|
|
|
CALL cp_calc_energy_force(f_env_id,1,ierr)
|
|
IF (ierr/=0) STOP "calc_force1"
|
|
CALL cp_get_energy(f_env_id,e0,ierr)
|
|
IF (ierr/=0) STOP "get_energy1"
|
|
WRITE(*,*) "e0= ",e0
|
|
CALL cp_get_force(f_env_id,f,SIZE(f),ierr)
|
|
IF (ierr/=0) STOP "get_force"
|
|
CALL cp_get_pos(f_env_id,pos,SIZE(pos),ierr)
|
|
IF (ierr/=0) STOP "get_pos"
|
|
f_m=0._dp
|
|
DO i=1,SIZE(f)
|
|
f_m=f_m+f(i)**2
|
|
END DO
|
|
f_m=SQRT(f_m/natom)
|
|
WRITE(*,*)"f_m",f_m
|
|
WRITE(*,*)
|
|
err_max=0._dp
|
|
DO i=1,SIZE(pos)
|
|
p0=pos(i)
|
|
pos(i)=p0+small
|
|
CALL cp_set_pos(f_env_id,pos,SIZE(pos),ierr)
|
|
IF (ierr/=0) STOP "set_pos1"
|
|
CALL cp_calc_energy_force(f_env_id,1,ierr)
|
|
IF (ierr/=0) STOP "calc_energy_force2"
|
|
CALL cp_get_energy(f_env_id,ep,ierr)
|
|
IF (ierr/=0) STOP "get_energy2"
|
|
pos(i)=p0-small
|
|
CALL cp_calc_energy(f_env_id,pos,SIZE(pos),em,ierr)
|
|
IF (ierr/=0) STOP "calc_energy"
|
|
pos(i)=p0
|
|
err=(em-ep)/(2._dp*small)-f(i)
|
|
WRITE(*,*) "err= ",err
|
|
if (abs(err)>abs(err_max)) err_max=err
|
|
WRITE (*,*) "f= ",f(i)," f_fdiff= ",(em-ep)/(2._dp*small),&
|
|
" err= ",err," err_m= ",err/f_m*100._dp
|
|
END DO
|
|
|
|
WRITE (*,*)
|
|
WRITE (*,*) "err_max ",err_max," err_max_m ",err_max/f_m*100._dp
|
|
!FM CALL cp_destroy_force_env(f_env_id,ierr)
|
|
!FM IF (ierr/=0) STOP "destroy_force_env"
|
|
|
|
CALL cp_finalize_cp2k(1,ierr)
|
|
IF (ierr/=0) STOP "finalize_cp2k"
|
|
|
|
END PROGRAM f_diff_low
|