Three pigeons with one bean: Bug fix for PBC in coorect bondend electrostatics both

energy/forces and virial.. Major rewriting of the routine and heavy cleaning.. Bug
fix in stress tensor for shell models (was missing the contribution from shell_i//core_i)..
Added two more regtests to check that PBC will be preserved since now on..
Resetted several regtests due to numerics (and to the bug fix for the shell model)


svn-origin-rev: 6124
This commit is contained in:
Teodoro Laino 2007-06-05 22:20:53 +00:00
parent 8028fe790d
commit 9b09991d00
15 changed files with 733 additions and 276 deletions

View file

@ -378,16 +378,16 @@ SUBROUTINE fist_force_control ( fist_env, virial, para_env, debug , force_env_s
IF(shell_present) THEN
CALL bonded_correct_gaussian ( atomic_kind_set, local_particles, &
particle_set, ewald_env, &
thermo % e_bonded, pv_bc ,&
shell_particle_set=shell_particle_set,&
core_particle_set=core_particle_set,&
error=error)
particle_set, ewald_env, &
thermo % e_bonded, pv_bc ,&
shell_particle_set=shell_particle_set,&
core_particle_set=core_particle_set,&
cell=cell, error=error)
ELSE
CALL bonded_correct_gaussian ( atomic_kind_set, local_particles, &
particle_set, ewald_env, &
thermo % e_bonded, pv_bc ,error=error)
particle_set, ewald_env, thermo % e_bonded,&
pv_bc=pv_bc , cell=cell, error=error)
END IF
IF ( first_time ) THEN
IF (iw > 0 ) CALL ewald_print ( iw , cell, vg_coulomb, &

View file

@ -2,6 +2,7 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 - 2007 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/fist_nonbond_force [1.0] *
!!
!! NAME
@ -26,7 +27,8 @@ MODULE fist_nonbond_force
USE atomic_kind_types, ONLY: atomic_kind_type,&
get_atomic_kind,&
get_atomic_kind_set
USE cell_types, ONLY: cell_type
USE cell_types, ONLY: cell_type,&
pbc
USE distribution_1d_types, ONLY: distribution_1d_type
USE ewald_environment_types, ONLY: ewald_env_get,&
ewald_environment_type
@ -61,21 +63,31 @@ MODULE fist_nonbond_force
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'fist_nonbond_force'
PUBLIC :: force_nonbond,&
bonded_correct_gaussian
bonded_correct_gaussian
!!***
!******************************************************************************
CONTAINS
!******************************************************************************
!!****** fist_nonbond_force/force_nonbond [1.0] *
!!
!! NAME
!! force_nonbond
!!
!! FUNCTION
!! Calculates the force and the potential of the minimum image, and
!! the pressure tensor
!!
!! AUTHOR
!!
!!
!! MODIFICATION HISTORY
!!
!!
!! SOURCE
!******************************************************************************
SUBROUTINE force_nonbond ( fist_nonbond_env, particle_set, cell, &
pot_nonbond, f_nonbond, ptens_nonbond, &
fshell_nonbond, fcore_nonbond, error )
! Calculates the force and the potential of the minimum image, and
! the pressure tensor
!
TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
TYPE(particle_type), DIMENSION(:), &
INTENT(IN) :: particle_set
@ -90,8 +102,8 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'force_nonbond', &
routineP = moduleN//':'//routineN
INTEGER :: atom_a, atom_b, handle, &
i, ikind, ilist, ipair, jkind, &
INTEGER :: atom_a, atom_b, handle, i, &
ikind, ilist, ipair, jkind, &
nkinds, npairs, shell_a, &
shell_b, shell_type
INTEGER, DIMENSION(:, :), POINTER :: list
@ -429,7 +441,7 @@ CONTAINS
pv_com(3,1) = pv_com(3,1) + rab(3) * fr(1)
pv_com(3,2) = pv_com(3,2) + rab(3) * fr(2)
pv_com(3,3) = pv_com(3,3) + rab(3) * fr(3)
!
END IF ! shell_type
END IF ! do_shell
ptens11 = ptens11 + pv_com(1,1)
@ -458,11 +470,26 @@ CONTAINS
CALL timestop ( handle )
END SUBROUTINE force_nonbond
!******************************************************************************
!!****** fist_nonbond_force/bonded_correct_gaussian [1.0] *
!!
!! NAME
!! bonded_correct_gaussian
!!
!! FUNCTION
!! corrects electrostatics for bonded terms
!!
!! AUTHOR
!!
!!
!! MODIFICATION HISTORY
!! Splitted routines to clean and to fix a bug with the tensor whose
!! original definition was not correct for PBC.. [Teodoro Laino -06/2007]
!!
!! SOURCE
!******************************************************************************
SUBROUTINE bonded_correct_gaussian ( atomic_kind_set, local_particles, &
particle_set, ewald_env, v_bonded_corr, pv_bc, shell_particle_set, &
core_particle_set, error )
core_particle_set, cell, error )
TYPE(atomic_kind_type), POINTER :: atomic_kind_set( : )
@ -474,12 +501,10 @@ CONTAINS
INTENT(OUT) :: pv_bc
TYPE(particle_type), OPTIONAL, POINTER :: shell_particle_set( : ), &
core_particle_set( : )
TYPE(cell_type), POINTER :: cell
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(LEN=*), PARAMETER :: routineN = 'bonded_correct_gaussian'
REAL(KIND=dp), PARAMETER :: ac1 = 0.254829592_dp, ac2 = -0.284496736_dp, &
ac3 = 1.421413741_dp, ac4 = -1.453152027_dp, ac5 = 1.061405429_dp, &
pc = 0.3275911_dp
INTEGER :: group, handle, i, iatom, &
iexl, ikind, j, &
@ -487,294 +512,101 @@ CONTAINS
shell_i_index, shell_j_index
LOGICAL :: do_shell, i_is_shell, &
j_is_shell, shell_adiabatic
REAL(KIND=dp) :: alpha, arg, const, dij, e_arg_arg, efac, eps0, errf, &
ffac, flops, fscalar, idij, ifourpieps0, qci, qcj, qi, qj, qsi, qsj, &
rijsq, tc, tc2
REAL(KIND=dp), DIMENSION(3) :: fij_com, rci, rcj, ri, rij, &
rj, rsi, rsj
TYPE(atomic_kind_type), POINTER :: atomic_kind
TYPE(shell_kind_type), POINTER :: shell_i, shell_j
! parameters for the numerical erf in Abramowitz and Stegun Eq. 7.1.26 pg 299
!------------------------------------------------------------------------------
REAL(KIND=dp) :: alpha, const, eps0, flops, &
ifourpieps0, qci, qcj, qsi, &
qsj
REAL(KIND=dp), DIMENSION(3) :: fij_com, rci, rcj, rsi, rsj
CALL timeset ( routineN, 'I', 'Mflops', handle )
flops = 0.0_dp
NULLIFY(atomic_kind, shell_i, shell_j)
! defining the constants
pv_bc = 0.0_dp
! Initializing values
pv_bc = 0.0_dp
v_bonded_corr = 0.0_dp
! Defining the constants
CALL ewald_env_get ( ewald_env, eps0 = eps0, alpha = alpha, group = group ,error=error)
const = 2.0_dp * alpha / SQRT( pi )
ifourpieps0 = 1.0_dp / ( 4.0_dp * pi * eps0 )
CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, shell_adiabatic=shell_adiabatic)
CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,&
shell_adiabatic=shell_adiabatic)
nkind = SIZE ( atomic_kind_set )
DO ikind = 1, nkind
natoms_per_kind = local_particles % n_el ( ikind )
DO iatom = 1, natoms_per_kind
i = local_particles % list ( ikind ) % array ( iatom )
atomic_kind => particle_set ( i ) % atomic_kind
ri ( : ) = particle_set ( i ) % r ( : )
CALL get_atomic_kind ( atomic_kind,qeff = qi,&
shell_active=i_is_shell,&
shell=shell_i )
IF(i_is_shell) THEN
shell_i_index = particle_set(i)%shell_index
rsi = shell_particle_set(shell_i_index)%r
rci = core_particle_set(shell_i_index)%r
qci = shell_i%charge_core
qsi = shell_i%charge_shell
qi = qci
ELSE
qsi = 0.0_dp
qci = qi
rci = ri
END IF
CALL get_atom_info(i, particle_set, shell_particle_set, core_particle_set,&
shell_i_index, i_is_shell, rsi, rci, qci, qsi)
DO iexl = 1, SIZE(particle_set(i)%list_exclude_ei)
j = particle_set(i)%list_exclude_ei(iexl)
IF ( j==i ) CYCLE
fij_com = 0.0_dp
atomic_kind => particle_set ( j ) % atomic_kind
rj ( : ) = particle_set ( j ) % r ( : )
CALL get_atomic_kind ( atomic_kind,qeff = qj,&
shell_active=j_is_shell, &
shell=shell_j )
IF(j_is_shell) THEN
shell_j_index = particle_set(j)%shell_index
rsj = shell_particle_set(shell_j_index)%r
rcj = core_particle_set(shell_j_index)%r
qcj = shell_j%charge_core
qsj = shell_j%charge_shell
qj = qcj
ELSE
qsj = 0.0_dp
qcj = qj
rcj = rj
END IF
j = particle_set(i)%list_exclude_ei(iexl)
IF ( j>=i ) CYCLE
CALL get_atom_info(j, particle_set, shell_particle_set, core_particle_set,&
shell_j_index, j_is_shell, rsj, rcj, qcj, qsj)
do_shell = .FALSE.
IF (i_is_shell.OR.j_is_shell) do_shell = .TRUE.
! shell-model: exclude core-core, core-ion, or ion-core interactions
! if not shell-model: exclude ion-ion interactions
rij(1) = rci(1) - rcj(1)
rij(2) = rci(2) - rcj(2)
rij(3) = rci(3) - rcj(3)
rijsq = rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3)
idij = 1.0_dp / SQRT ( rijsq )
dij = rijsq * idij
arg = alpha * dij
e_arg_arg = EXP ( -arg ** 2 )
tc = 1.0_dp / ( 1.0_dp + pc * arg )
! defining errf=1-erfc
errf = 1.0_dp &
- ((((ac5*tc+ac4)*tc+ac3)*tc+ac2)*tc+ac1) * tc * e_arg_arg
! getting the potential
v_bonded_corr = v_bonded_corr - qci*qcj*idij*errf
! subtracting the force from the total force
fscalar = ifourpieps0 * qci * qcj * idij ** 2 &
* ( idij * errf - const * e_arg_arg )
!
! This term excludes always ION-ION interactions
IF(i_is_shell) THEN
core_particle_set ( shell_i_index ) % f(1) = &
core_particle_set ( shell_i_index ) % f(1) - fscalar*rij(1)
core_particle_set ( shell_i_index ) % f(2) = &
core_particle_set ( shell_i_index ) % f(2) - fscalar*rij(2)
core_particle_set ( shell_i_index ) % f(3) = &
core_particle_set ( shell_i_index ) % f(3) - fscalar*rij(3)
CALL bonded_correct_gaussian_low(rci, rcj, cell, v_bonded_corr, &
core_particle_set, core_particle_set, &
shell_i_index, shell_j_index, .TRUE., &
alpha, qci, qcj, ifourpieps0, const, pv_bc)
ELSE
particle_set ( i ) % f(1) = particle_set ( i ) % f(1) - fscalar*rij(1)
particle_set ( i ) % f(2) = particle_set ( i ) % f(2) - fscalar*rij(2)
particle_set ( i ) % f(3) = particle_set ( i ) % f(3) - fscalar*rij(3)
CALL bonded_correct_gaussian_low(rci, rcj, cell, v_bonded_corr,&
particle_set, particle_set, i, j, .TRUE.,&
alpha, qci, qcj, ifourpieps0, const, pv_bc)
END IF
fij_com = fscalar*rij
! shell-model: exclude shell-shell interactions
! Shell-Model
IF (do_shell) THEN
! shell-model: exclude shell_i-shell_j interactions
IF(i_is_shell .AND. j_is_shell) THEN
rij(1) = rsi(1) - rsj(1)
rij(2) = rsi(2) - rsj(2)
rij(3) = rsi(3) - rsj(3)
rijsq = rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3)
idij = 1.0_dp / SQRT ( rijsq )
dij = rijsq * idij
arg = alpha * dij
e_arg_arg = EXP ( -arg ** 2 )
tc = 1.0_dp / ( 1.0_dp + pc * arg )
! defining errf=1-erfc
errf = 1.0_dp &
- ((((ac5*tc+ac4)*tc+ac3)*tc+ac2)*tc+ac1) * tc * e_arg_arg
! getting the potential
v_bonded_corr = v_bonded_corr - qsi*qsj*idij*errf
! subtracting the force from the total force
fscalar = ifourpieps0 * qsi * qsj * idij ** 2 &
* ( idij * errf - const * e_arg_arg )
shell_particle_set ( shell_i_index) % f(1) = &
shell_particle_set ( shell_i_index ) % f(1) - fscalar*rij(1)
shell_particle_set ( shell_i_index ) % f(2) = &
shell_particle_set ( shell_i_index ) % f(2) - fscalar*rij(2)
shell_particle_set ( shell_i_index ) % f(3) = &
shell_particle_set ( shell_i_index ) % f(3) - fscalar*rij(3)
IF(shell_adiabatic) THEN
fij_com = fij_com + fscalar*rij
END IF
CALL bonded_correct_gaussian_low(rsi, rsj, cell, v_bonded_corr,&
shell_particle_set, shell_particle_set, shell_i_index, shell_j_index,&
shell_adiabatic, alpha, qsi, qsj, ifourpieps0, const, pv_bc)
END IF
! shell-model: exclude shell_i-core interactions
! shell-model: exclude shell_i-core_j interactions
IF(i_is_shell ) THEN
rij(1) = rsi(1) - rcj(1)
rij(2) = rsi(2) - rcj(2)
rij(3) = rsi(3) - rcj(3)
rijsq = rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3)
idij = 1.0_dp / SQRT ( rijsq )
dij = rijsq * idij
arg = alpha * dij
e_arg_arg = EXP ( -arg ** 2 )
tc = 1.0_dp / ( 1.0_dp + pc * arg )
! defining errf=1-erfc
errf = 1.0_dp &
- ((((ac5*tc+ac4)*tc+ac3)*tc+ac2)*tc+ac1) * tc * e_arg_arg
! getting the potential
v_bonded_corr = v_bonded_corr - qsi*qcj*idij*errf
! subtracting the force from the total force
fscalar = ifourpieps0 * qsi * qcj * idij ** 2 &
* ( idij * errf - const * e_arg_arg )
shell_particle_set ( shell_i_index) % f(1) = &
shell_particle_set ( shell_i_index ) % f(1) - fscalar*rij(1)
shell_particle_set ( shell_i_index ) % f(2) = &
shell_particle_set ( shell_i_index ) % f(2) - fscalar*rij(2)
shell_particle_set ( shell_i_index ) % f(3) = &
shell_particle_set ( shell_i_index ) % f(3) - fscalar*rij(3)
IF(shell_adiabatic) THEN
fij_com = fij_com + fscalar*rij
END IF
END IF
IF(j_is_shell ) THEN
rij(1) = rci(1) - rsj(1)
rij(2) = rci(2) - rsj(2)
rij(3) = rci(3) - rsj(3)
rijsq = rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3)
idij = 1.0_dp / SQRT ( rijsq )
dij = rijsq * idij
arg = alpha * dij
e_arg_arg = EXP ( -arg ** 2 )
tc = 1.0_dp / ( 1.0_dp + pc * arg )
! defining errf=1-erfc
errf = 1.0_dp &
- ((((ac5*tc+ac4)*tc+ac3)*tc+ac2)*tc+ac1) * tc * e_arg_arg
! getting the potential
v_bonded_corr = v_bonded_corr - qci*qsj*idij*errf
! subtracting the force from the total force
fscalar = ifourpieps0 * qci * qsj * idij ** 2 &
* ( idij * errf - const * e_arg_arg )
IF(i_is_shell) THEN
core_particle_set ( shell_i_index) % f(1) = &
core_particle_set ( shell_i_index ) % f(1) - fscalar*rij(1)
core_particle_set ( i) % f(2) = &
core_particle_set ( shell_i_index ) % f(2) - fscalar*rij(2)
core_particle_set (i) % f(3) = &
core_particle_set ( shell_i_index ) % f(3) - fscalar*rij(3)
IF (j_is_shell) THEN
CALL bonded_correct_gaussian_low(rsi, rcj, cell, v_bonded_corr,&
shell_particle_set, core_particle_set, shell_i_index, shell_j_index,&
shell_adiabatic, alpha, qsi, qcj, ifourpieps0, const, pv_bc)
ELSE
particle_set ( i) % f(1) = &
particle_set ( i ) % f(1) - fscalar*rij(1)
particle_set ( i) % f(2) = &
particle_set ( i ) % f(2) - fscalar*rij(2)
particle_set (i) % f(3) = &
particle_set ( i ) % f(3) - fscalar*rij(3)
CALL bonded_correct_gaussian_low(rsi, rcj, cell, v_bonded_corr,&
shell_particle_set, particle_set, shell_i_index, j, shell_adiabatic,&
alpha, qsi, qcj, ifourpieps0, const, pv_bc)
END IF
IF(shell_adiabatic) THEN
fij_com = fij_com + fscalar*rij
END IF
! shell-model: exclude shell_j-core_i interactions
IF(j_is_shell ) THEN
IF (i_is_shell) THEN
CALL bonded_correct_gaussian_low(rci, rsj, cell, v_bonded_corr,&
core_particle_set, shell_particle_set, shell_i_index, shell_j_index,&
shell_adiabatic, alpha, qci, qsj, ifourpieps0, const, pv_bc)
ELSE
CALL bonded_correct_gaussian_low(rci, rsj, cell, v_bonded_corr,&
particle_set, shell_particle_set, i, shell_j_index, shell_adiabatic,&
alpha, qci, qsj, ifourpieps0, const, pv_bc)
END IF
END IF
END IF
pv_bc(1,1) = pv_bc(1,1) - fij_com(1) * particle_set ( i ) % r(1)
pv_bc(1,2) = pv_bc(1,2) - fij_com(1) * particle_set ( i ) % r(2)
pv_bc(1,3) = pv_bc(1,3) - fij_com(1) * particle_set ( i ) % r(3)
pv_bc(2,1) = pv_bc(2,1) - fij_com(2) * particle_set ( i ) % r(1)
pv_bc(2,2) = pv_bc(2,2) - fij_com(2) * particle_set ( i ) % r(2)
pv_bc(2,3) = pv_bc(2,3) - fij_com(2) * particle_set ( i ) % r(3)
pv_bc(3,1) = pv_bc(3,1) - fij_com(3) * particle_set ( i ) % r(1)
pv_bc(3,2) = pv_bc(3,2) - fij_com(3) * particle_set ( i ) % r(2)
pv_bc(3,3) = pv_bc(3,3) - fij_com(3) * particle_set ( i ) % r(3)
flops = flops + 62.0_dp
ENDDO
! Exclusion of shell-core interaction
! Always Exclude shell_i-core_i interaction
IF(i_is_shell ) THEN
rij(1) = rci(1) - rsi(1)
rij(2) = rci(2) - rsi(2)
rij(3) = rci(3) - rsi(3)
rijsq = rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3)
dij = SQRT(rijsq)
! Two possible limiting cases according the value of dij
arg = alpha * dij
! and this is a magic number.. it is related to the order expansion
! and to the value of the polynomial coefficients
IF(arg > 0.355_dp) THEN
idij = 1.0_dp / dij
e_arg_arg = EXP ( -arg ** 2 )
tc = 1.0_dp / ( 1.0_dp + pc * arg )
! defining errf=1-erfc
errf = 1.0_dp &
- ((((ac5*tc+ac4)*tc+ac3)*tc+ac2)*tc+ac1) * tc * e_arg_arg
efac = idij*errf
ffac = idij**2 * (idij*errf - const*e_arg_arg)
ELSE
tc = arg*arg
tc2 = tc*tc
efac = 1.0_dp - tc/3.0_dp + tc2/10.0_dp - tc*tc2/42.0_dp + tc2*tc2/216.0_dp
efac = efac * const
ffac = - 2.0_dp/3.0_dp + 2.0_dp/5.0_dp*tc - 1.0_dp/7.0_dp*tc2 + 1.0_dp/27.0_dp*tc2*tc
ffac = const*alpha**2*ffac
ENDIF
! getting the potential
! th factor 2 comes to compensate for the factor 0.5
! introduced later to avoid double counting
v_bonded_corr = v_bonded_corr - 2.0_dp * qci*qsi*efac
! subtracting the force from the total force
fscalar = ifourpieps0 * qci * qsi * ffac
core_particle_set ( shell_i_index) % f(1) = &
core_particle_set ( shell_i_index ) % f(1) - fscalar*rij(1)
core_particle_set ( shell_i_index) % f(2) = &
core_particle_set ( shell_i_index ) % f(2) - fscalar*rij(2)
core_particle_set (shell_i_index) % f(3) = &
core_particle_set ( shell_i_index ) % f(3) - fscalar*rij(3)
shell_particle_set (shell_i_index) % f(1) = &
shell_particle_set ( shell_i_index ) % f(1) + fscalar*rij(1)
shell_particle_set ( shell_i_index) % f(2) = &
shell_particle_set ( shell_i_index ) % f(2) + fscalar*rij(2)
shell_particle_set (shell_i_index) % f(3) = &
shell_particle_set ( shell_i_index ) % f(3) + fscalar*rij(3)
CALL bonded_correct_gaussian_low_sh(rci, rsi, cell, v_bonded_corr,&
core_particle_set, shell_particle_set, shell_i_index, shell_adiabatic,&
alpha, qci, qsi, ifourpieps0, const, pv_bc)
END IF
END DO
ENDDO
! the factor of 1/2 comes from double counting in the exclusion list
v_bonded_corr = v_bonded_corr * ifourpieps0 * 0.5_dp
v_bonded_corr = v_bonded_corr * ifourpieps0
CALL mp_sum ( v_bonded_corr, group )
flops = flops * 1.0E-6_dp
@ -782,6 +614,237 @@ CONTAINS
END SUBROUTINE bonded_correct_gaussian
!******************************************************************************
!!****** fist_nonbond_force/bonded_correct_gaussian_low [1.0] *
!!
!! NAME
!! bonded_correct_gaussian_low
!!
!! FUNCTION
!!
!! AUTHOR
!! Teodoro Laino
!!
!! MODIFICATION HISTORY
!! Splitted routines to clean and to fix a bug with the tensor whose
!! original definition was not correct for PBC..
!!
!! SOURCE
!******************************************************************************
SUBROUTINE bonded_correct_gaussian_low(r1, r2, cell, v_bonded_corr,&
particle_set1, particle_set2, i, j, shell_adiabatic, alpha, q1, q2, &
ifourpieps0, const, pv_bc)
REAL(KIND=dp), DIMENSION(3) :: r1, r2
TYPE(cell_type), POINTER :: cell
REAL(KIND=dp), INTENT(INOUT) :: v_bonded_corr
TYPE(particle_type), POINTER :: particle_set1(:), &
particle_set2(:)
INTEGER, INTENT(IN) :: i, j
LOGICAL, INTENT(IN) :: shell_adiabatic
REAL(KIND=dp), INTENT(IN) :: alpha, q1, q2, ifourpieps0, &
const
REAL(KIND=dp), INTENT(INOUT) :: pv_bc(3,3)
CHARACTER(LEN=*), PARAMETER :: routineN = 'bonded_correct_gaussian_low'
REAL(KIND=dp), PARAMETER :: ac1 = 0.254829592_dp, ac2 = -0.284496736_dp, &
ac3 = 1.421413741_dp, ac4 = -1.453152027_dp, ac5 = 1.061405429_dp, &
pc = 0.3275911_dp
REAL(KIND=dp) :: arg, dij, e_arg_arg, errf, &
fscalar, idij, rijsq, tc
REAL(KIND=dp), DIMENSION(3) :: fij_com, rij
rij = r1 - r2
rij = pbc(rij, cell)
rijsq = rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3)
idij = 1.0_dp / SQRT ( rijsq )
dij = rijsq * idij
arg = alpha * dij
e_arg_arg = EXP ( -arg ** 2 )
tc = 1.0_dp / ( 1.0_dp + pc * arg )
! Defining errf=1-erfc
errf = 1.0_dp &
- ((((ac5*tc+ac4)*tc+ac3)*tc+ac2)*tc+ac1) * tc * e_arg_arg
! Getting the potential
v_bonded_corr = v_bonded_corr - q1*q2*idij*errf
! Subtracting the force from the total force
fscalar = ifourpieps0 * q1 * q2 * idij**2 * (idij*errf-const*e_arg_arg)
particle_set1(i)%f(1) = particle_set1(i)%f(1) - fscalar*rij(1)
particle_set1(i)%f(2) = particle_set1(i)%f(2) - fscalar*rij(2)
particle_set1(i)%f(3) = particle_set1(i)%f(3) - fscalar*rij(3)
particle_set2(j)%f(1) = particle_set2(j)%f(1) + fscalar*rij(1)
particle_set2(j)%f(2) = particle_set2(j)%f(2) + fscalar*rij(2)
particle_set2(j)%f(3) = particle_set2(j)%f(3) + fscalar*rij(3)
IF(shell_adiabatic) THEN
fij_com = fscalar*rij
pv_bc(1,1) = pv_bc(1,1) - fij_com(1) * rij(1)
pv_bc(1,2) = pv_bc(1,2) - fij_com(1) * rij(2)
pv_bc(1,3) = pv_bc(1,3) - fij_com(1) * rij(3)
pv_bc(2,1) = pv_bc(2,1) - fij_com(2) * rij(1)
pv_bc(2,2) = pv_bc(2,2) - fij_com(2) * rij(2)
pv_bc(2,3) = pv_bc(2,3) - fij_com(2) * rij(3)
pv_bc(3,1) = pv_bc(3,1) - fij_com(3) * rij(1)
pv_bc(3,2) = pv_bc(3,2) - fij_com(3) * rij(2)
pv_bc(3,3) = pv_bc(3,3) - fij_com(3) * rij(3)
END IF
END SUBROUTINE bonded_correct_gaussian_low
!!****** fist_nonbond_force/bonded_correct_gaussian_low_sh [1.0] *
!!
!! NAME
!! bonded_correct_gaussian_low_sh
!!
!! FUNCTION
!! specific for shell models cleans the interaction core-shell on the same
!! atom
!!
!! AUTHOR
!! Teodoro Laino
!!
!! MODIFICATION HISTORY
!! Splitted routines to clean and to fix a bug with the tensor whose
!! original definition was not correct for PBC..
!!
!! SOURCE
!******************************************************************************
SUBROUTINE bonded_correct_gaussian_low_sh(r1, r2, cell, v_bonded_corr, &
core_particle_set, shell_particle_set, i, shell_adiabatic, alpha, q1, q2,&
ifourpieps0, const, pv_bc)
REAL(KIND=dp), DIMENSION(3) :: r1, r2
TYPE(cell_type), POINTER :: cell
REAL(KIND=dp), INTENT(INOUT) :: v_bonded_corr
TYPE(particle_type), POINTER :: core_particle_set(:), &
shell_particle_set(:)
INTEGER, INTENT(IN) :: i
LOGICAL, INTENT(IN) :: shell_adiabatic
REAL(KIND=dp), INTENT(IN) :: alpha, q1, q2, ifourpieps0, &
const
REAL(KIND=dp), INTENT(INOUT) :: pv_bc(3,3)
CHARACTER(LEN=*), PARAMETER :: routineN = 'bonded_correct_gaussian_low_sh'
REAL(KIND=dp), PARAMETER :: ac1 = 0.254829592_dp, ac2 = -0.284496736_dp, &
ac3 = 1.421413741_dp, ac4 = -1.453152027_dp, ac5 = 1.061405429_dp, &
pc = 0.3275911_dp
REAL(KIND=dp) :: arg, dij, e_arg_arg, efac, &
errf, ffac, fscalar, idij, &
rijsq, tc, tc2
REAL(KIND=dp), DIMENSION(3) :: fij_com, rij
rij = r1 - r2
rij = pbc(rij, cell)
rijsq = rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3)
dij = SQRT(rijsq)
! Two possible limiting cases according the value of dij
arg = alpha * dij
! and this is a magic number.. it is related to the order expansion
! and to the value of the polynomial coefficients
IF(arg > 0.355_dp) THEN
idij = 1.0_dp / dij
e_arg_arg = EXP ( -arg ** 2 )
tc = 1.0_dp / ( 1.0_dp + pc * arg )
! defining errf=1-erfc
errf = 1.0_dp &
- ((((ac5*tc+ac4)*tc+ac3)*tc+ac2)*tc+ac1) * tc * e_arg_arg
efac = idij*errf
ffac = idij**2 * (idij*errf - const*e_arg_arg)
ELSE
tc = arg*arg
tc2 = tc*tc
efac = 1.0_dp - tc/3.0_dp + tc2/10.0_dp - tc*tc2/42.0_dp + tc2*tc2/216.0_dp
efac = efac * const
ffac = - 2.0_dp/3.0_dp + 2.0_dp/5.0_dp*tc - 1.0_dp/7.0_dp*tc2 + 1.0_dp/27.0_dp*tc2*tc
ffac = const*alpha**2*ffac
ENDIF
! getting the potential
v_bonded_corr = v_bonded_corr - q1*q2*efac
! subtracting the force from the total force
fscalar = ifourpieps0 * q1 * q2 * ffac
core_particle_set(i)%f(1)=core_particle_set(i)%f(1)-fscalar*rij(1)
core_particle_set(i)%f(2)=core_particle_set(i)%f(2)-fscalar*rij(2)
core_particle_set(i)%f(3)=core_particle_set(i)%f(3)-fscalar*rij(3)
shell_particle_set(i)%f(1)=shell_particle_set(i)%f(1)+fscalar*rij(1)
shell_particle_set(i)%f(2)=shell_particle_set(i)%f(2)+fscalar*rij(2)
shell_particle_set(i)%f(3)=shell_particle_set(i)%f(3)+fscalar*rij(3)
IF(shell_adiabatic) THEN
fij_com = fscalar*rij
pv_bc(1,1) = pv_bc(1,1) - fij_com(1) * rij(1)
pv_bc(1,2) = pv_bc(1,2) - fij_com(1) * rij(2)
pv_bc(1,3) = pv_bc(1,3) - fij_com(1) * rij(3)
pv_bc(2,1) = pv_bc(2,1) - fij_com(2) * rij(1)
pv_bc(2,2) = pv_bc(2,2) - fij_com(2) * rij(2)
pv_bc(2,3) = pv_bc(2,3) - fij_com(2) * rij(3)
pv_bc(3,1) = pv_bc(3,1) - fij_com(3) * rij(1)
pv_bc(3,2) = pv_bc(3,2) - fij_com(3) * rij(2)
pv_bc(3,3) = pv_bc(3,3) - fij_com(3) * rij(3)
END IF
END SUBROUTINE bonded_correct_gaussian_low_sh
!!****** fist_nonbond_force/get_atom_info [1.0] *
!!
!! NAME
!! get_atom_info
!!
!! FUNCTION
!! gives back information on atoms.. specific for bonded correct gaussian
!!
!! AUTHOR
!! Teodoro Laino
!!
!! MODIFICATION HISTORY
!! Splitted routines to clean and to fix a bug with the tensor whose
!! original definition was not correct for PBC..
!!
!! SOURCE
!******************************************************************************
SUBROUTINE get_atom_info(i, particle_set, shell_particle_set, core_particle_set,&
shell_i_index, i_is_shell, rsi, rci, qci, qsi)
INTEGER, INTENT(IN) :: i
TYPE(particle_type), POINTER :: particle_set(:)
TYPE(particle_type), OPTIONAL, POINTER :: shell_particle_set(:), &
core_particle_set(:)
INTEGER, INTENT(OUT) :: shell_i_index
LOGICAL, INTENT(OUT) :: i_is_shell
REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: rsi, rci
REAL(KIND=dp), INTENT(OUT) :: qci, qsi
CHARACTER(LEN=*), PARAMETER :: routineN = 'get_atom_info'
REAL(KIND=dp) :: qi
REAL(KIND=dp), DIMENSION(3) :: ri
TYPE(atomic_kind_type), POINTER :: atomic_kind
TYPE(shell_kind_type), POINTER :: shell_i
NULLIFY(atomic_kind, shell_i)
atomic_kind => particle_set ( i ) % atomic_kind
ri (:) = particle_set ( i ) % r (:)
CALL get_atomic_kind ( atomic_kind,qeff = qi,&
shell_active=i_is_shell,&
shell=shell_i )
IF(i_is_shell) THEN
shell_i_index = particle_set(i)%shell_index
rsi = shell_particle_set(shell_i_index)%r
rci = core_particle_set(shell_i_index)%r
qci = shell_i%charge_core
qsi = shell_i%charge_shell
ELSE
qsi = 0.0_dp
qci = qi
rci = ri
END IF
END SUBROUTINE get_atom_info
END MODULE fist_nonbond_force

View file

@ -5,3 +5,10 @@
argon_GENPOT_fcc-ql.inp
# reset for smearing
argon_GENPOT_fcc-ql.inp
# faster bonded correct gaussian and bug fix for stress tensor in shell model
water_3_dist_2.inp
# faster bonded correct gaussian and bug fix for stress tensor in shell model
water_3_dist_3.inp
water_3_dist_inter2.inp
# faster bonded correct gaussian and bug fix for stress tensor in shell model
water_3_dist_startC.inp

View file

@ -31,3 +31,6 @@ water_3_dist_g.inp 2
water_3_dist_centers.inp 2
water_3_dist_centers2.inp 2
water_3_dist_centers3.inp 2
# more on pbc
cubane_15.inp 7
cubane_0.inp 7

View file

@ -13,3 +13,8 @@ uo2_shell_npt300.inp
C_tersoff.inp
#numerics only for online regtest
water_3_dist_centers2.inp
# faster bonded correct gaussian and bug fix for stress tensor in shell model
uo2_shell_npt300.inp
uo2_shell_npti_b.inp
uo2_shell_npti.inp
water_3_dist_g.inp

View file

@ -0,0 +1,66 @@
&FORCE_EVAL
METHOD FIST
&MM
&FORCEFIELD
! PARM91
! VDW_SCALE14 0.5
! EI_SCALE14 0.5
! PARM94 and Following
VDW_SCALE14 0.5
EI_SCALE14 0.83333333333333333333
parm_file_name ../sample_pot/cubane.pot
parmtype CHM
scale_cutoff .false.
&END FORCEFIELD
&POISSON
&EWALD
EWALD_TYPE spme
ALPHA .36
GMAX 64
&END EWALD
&END POISSON
&END MM
&SUBSYS
&CELL
ABC 30.0 30.0 30.0
UNIT ANGSTROM
&END CELL
&TOPOLOGY
CONN_FILE_NAME ../sample_psf/cubane.psf
CONNECTIVITY PSF
COORD_FILE_NAME ../sample_pdb/cubane_0.pdb
COORDINATE PDB
&END TOPOLOGY
&PRINT
&TOPOLOGY_INFO
UTIL_INFO
&END
&END
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
FFTLIB FFTSG
PRINT_LEVEL DEBUG
PROJECT CUBANE
RUN_TYPE GEO_OPT
#RUN_TYPE MD
&END GLOBAL
&MOTION
&GEOOPT
MAX_ITER 50
&END
&MD
ENSEMBLE NPT_F
STEPS 5
TIMESTEP 0.5
&END
&PRINT
&TRAJECTORY
FILENAME __STD_OUT__
FORMAT ATOMIC
&END
&END
&END MOTION

View file

@ -0,0 +1,66 @@
&FORCE_EVAL
METHOD FIST
&MM
&FORCEFIELD
! PARM91
! VDW_SCALE14 0.5
! EI_SCALE14 0.5
! PARM94 and Following
VDW_SCALE14 0.5
EI_SCALE14 0.83333333333333333333
parm_file_name ../sample_pot/cubane.pot
parmtype CHM
scale_cutoff .false.
&END FORCEFIELD
&POISSON
&EWALD
EWALD_TYPE spme
ALPHA .36
GMAX 64
&END EWALD
&END POISSON
&END MM
&SUBSYS
&CELL
ABC 30.0 30.0 30.0
UNIT ANGSTROM
&END CELL
&TOPOLOGY
CONN_FILE_NAME ../sample_psf/cubane.psf
CONNECTIVITY PSF
COORD_FILE_NAME ../sample_pdb/cubane_15.pdb
COORDINATE PDB
&END TOPOLOGY
&PRINT
&TOPOLOGY_INFO
UTIL_INFO
&END
&END
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
FFTLIB FFTSG
PRINT_LEVEL DEBUG
PROJECT CUBANE
RUN_TYPE GEO_OPT
#RUN_TYPE MD
&END GLOBAL
&MOTION
&GEOOPT
MAX_ITER 50
&END
&MD
ENSEMBLE NPT_F
STEPS 5
TIMESTEP 0.5
&END
&PRINT
&TRAJECTORY
FILENAME __STD_OUT__
FORMAT ATOMIC
&END
&END
&END MOTION

View file

@ -273,3 +273,5 @@ water_3_rg3x3_ext.inp
C_tersoff.inp
# reduce the number of steps to keep runtime reasonable
C_tersoff.inp
# faster bonded correct gaussian and bug fix for stress tensor in shell model
water_3_dist.inp

View file

@ -0,0 +1,18 @@
CRYST1 30.000 30.000 30.000 90.00 90.00 90.00 P 1 1
ATOM 1 C1 MOL X 1 0.783 0.782 0.783 0.00 0.00 C
ATOM 2 H8 MOL X 1 1.408 1.407 1.408 0.00 0.00 H
ATOM 3 C2 MOL X 1 0.783 29.217 0.783 0.00 0.00 C
ATOM 4 H6 MOL X 1 1.408 28.602 1.408 0.00 0.00 H
ATOM 5 C3 MOL X 1 0.783 29.217 29.218 0.00 0.00 C
ATOM 6 H1 MOL X 1 1.408 28.592 28.593 0.00 0.00 H
ATOM 7 C4 MOL X 1 0.783 0.782 29.218 0.00 0.00 C
ATOM 8 H2 MOL X 1 1.408 1.407 28.593 0.00 0.00 H
ATOM 9 C5 MOL X 1 29.218 0.782 29.218 0.00 0.00 C
ATOM 10 H3 MOL X 1 28.583 1.407 28.593 0.00 0.00 H
ATOM 11 C6 MOL X 1 29.218 0.782 0.783 0.00 0.00 C
ATOM 12 H4 MOL X 1 28.583 1.407 1.408 0.00 0.00 H
ATOM 13 C7 MOL X 1 29.218 29.217 0.783 0.00 0.00 C
ATOM 14 H5 MOL X 1 28.583 28.592 1.408 0.00 0.00 H
ATOM 15 C8 MOL X 1 29.218 29.217 29.218 0.00 0.00 C
ATOM 16 H7 MOL X 1 28.583 28.592 28.593 0.00 0.00 H
END

View file

@ -0,0 +1,19 @@
REMARK MOL
ATOM 1 C1 MOL 1 15.783 15.782 15.783 C
ATOM 2 H8 MOL 1 16.408 16.407 16.408 H
ATOM 3 C2 MOL 1 15.783 14.217 15.783 C
ATOM 4 H6 MOL 1 16.408 13.602 16.408 H
ATOM 5 C3 MOL 1 15.783 14.217 14.218 C
ATOM 6 H1 MOL 1 16.408 13.592 13.593 H
ATOM 7 C4 MOL 1 15.783 15.782 14.218 C
ATOM 8 H2 MOL 1 16.408 16.407 13.593 H
ATOM 9 C5 MOL 1 14.218 15.782 14.218 C
ATOM 10 H3 MOL 1 13.583 16.407 13.593 H
ATOM 11 C6 MOL 1 14.218 15.782 15.783 C
ATOM 12 H4 MOL 1 13.583 16.407 16.408 H
ATOM 13 C7 MOL 1 14.218 14.217 15.783 C
ATOM 14 H5 MOL 1 13.583 13.592 16.408 H
ATOM 15 C8 MOL 1 14.218 14.217 14.218 C
ATOM 16 H7 MOL 1 13.583 13.592 13.593 H
TER
END

View file

@ -0,0 +1,81 @@
*>>>>>>> AMBER FF Converted into CHARMM FF style <<<<<<<
*>>>>>>> Generated on :: 20070605 124119.566 +0200 by :: teo <<<<<<<
*>>>>>>> pcihpc07 <<<<<<<
*>>>>>>> Leap Title :: MOL <<<<<<<
*>>>>>>> Send all comments related to the FFs conversion to <<<<<<<
*>>>>>>> teodoro.laino@gmail.com <<<<<<<
BONDS
!
!V(bond) = Kb(b - b0)**2
!
!Kb: kcal/mole/A**2
!b0: A
!
!atom type Kb b0
!
cy hc 334.500000000 1.094000000
cy cy 286.800000000 1.554000000
ANGLES
!
!V(angle) = Ktheta(Theta - Theta0)**2
!
!V(Urey-Bradley) = Kub(S - S0)**2
!
!Ktheta: kcal/mole/rad**2
!Theta0: degrees
!Kub: kcal/mole/A**2 (Urey-Bradley)
!S0: A
!
!atom types Ktheta Theta0 Kub S0
!
cy cy hc 44.800000000 115.140049232
cy cy cy 70.200000000 87.610037376
DIHEDRALS
!
!V(dihedral) = Kchi(1 + cos(n(chi) - delta))
!
!Kchi: kcal/mole
!n: multiplicity
!delta: degrees
!
!atom types Kchi n delta
!
cy cy cy hc 0.155555556 3 0.000000000
hc cy cy hc 0.155555556 3 0.000000000
cy cy cy cy 0.155555556 3 0.000000000
IMPROPER
!
!V(improper) = Kpsi(psi - psi0)**2
!
!Kpsi: kcal/mole/rad**2
!psi0: degrees
!note that the second column of numbers (0) is ignored
!
!atom types Kpsi psi0
!
NONBONDED
!
!V(Lennard-Jones) = Eps,i,j[(Rmin,i,j/ri,j)**12 - 2(Rmin,i,j/ri,j)**6]
!
!epsilon: kcal/mole, Eps,i,j = sqrt(eps,i * eps,j)
!Rmin/2: A, Rmin,i,j = Rmin/2,i + Rmin/2,j
!
!atom ignored epsilon Rmin/2 ignored eps,1-4 Rmin/2,1-4
!
cy 0.000000000 0.086000000 1.908000000
hc 0.000000000 0.015700000 1.487000000
END
!
! This Section can be cutted & pasted into the Fist input file..
!
CHARGES
cy -0.046700000
hc 0.046700000
END CHARGES

View file

@ -0,0 +1,114 @@
PSF EXT
1 !NTITLE
Conversion from AMBER PARMTOP ::MOL
16 !NATOM
1 MOL01 1 MOL C1 cy -0.046700 12.010 0
2 MOL01 1 MOL H8 hc 0.046700 1.008 0
3 MOL01 1 MOL C2 cy -0.046700 12.010 0
4 MOL01 1 MOL H6 hc 0.046700 1.008 0
5 MOL01 1 MOL C3 cy -0.046700 12.010 0
6 MOL01 1 MOL H1 hc 0.046700 1.008 0
7 MOL01 1 MOL C4 cy -0.046700 12.010 0
8 MOL01 1 MOL H2 hc 0.046700 1.008 0
9 MOL01 1 MOL C5 cy -0.046700 12.010 0
10 MOL01 1 MOL H3 hc 0.046700 1.008 0
11 MOL01 1 MOL C6 cy -0.046700 12.010 0
12 MOL01 1 MOL H4 hc 0.046700 1.008 0
13 MOL01 1 MOL C7 cy -0.046700 12.010 0
14 MOL01 1 MOL H5 hc 0.046700 1.008 0
15 MOL01 1 MOL C8 cy -0.046700 12.010 0
16 MOL01 1 MOL H7 hc 0.046700 1.008 0
20 !NBOND
1 2 3 4 5 6 7 8
9 10 11 12 13 14 15 16
1 3 1 7 1 11 3 5
3 13 5 7 5 15 7 9
9 11 9 15 11 13 13 15
48 !NTHETA
1 3 4 1 7 8 1 11 12
2 1 3 2 1 7 2 1 11
3 5 6 3 13 14 4 3 5
4 3 13 5 7 8 5 15 16
6 5 7 6 5 15 7 9 10
8 7 9 9 11 12 9 15 16
10 9 11 10 9 15 11 13 14
12 11 13 13 15 16 14 13 15
1 3 5 1 3 13 1 7 5
1 7 9 1 11 9 1 11 13
3 1 7 3 1 11 3 5 7
3 5 15 3 13 11 3 13 15
5 3 13 5 7 9 5 15 9
5 15 13 7 1 11 7 5 15
7 9 11 7 9 15 9 11 13
9 15 13 11 9 15 11 13 15
108 !NPHI
1 3 5 6 1 3 13 14
1 7 5 6 1 7 9 10
1 11 9 10 1 11 13 14
2 1 3 4 2 1 3 5
2 1 3 13 2 1 7 5
2 1 7 8 2 1 7 9
2 1 11 9 2 1 11 12
2 1 11 13 3 1 7 8
3 1 11 12 3 5 7 8
3 5 15 16 3 13 11 12
3 13 15 16 7 1 3 4
11 1 3 4 4 3 5 6
4 3 5 7 4 3 5 15
4 3 13 11 4 3 13 14
4 3 13 15 5 3 13 14
5 7 9 10 5 15 9 10
5 15 13 14 6 5 3 13
6 5 7 8 6 5 7 9
6 5 15 9 6 5 15 13
6 5 15 16 7 1 11 12
7 5 15 16 7 9 11 12
7 9 15 16 11 1 7 8
8 7 5 15 8 7 9 10
8 7 9 11 8 7 9 15
9 11 13 14 9 15 13 14
10 9 11 12 10 9 11 13
10 9 15 13 10 9 15 16
11 9 15 16 11 13 15 16
12 11 9 15 12 11 13 14
12 11 13 15 14 13 15 16
1 3 5 7 1 3 5 15
1 3 13 11 1 3 13 15
1 7 5 3 1 7 5 15
1 7 9 11 1 7 9 15
1 11 9 7 1 11 9 15
1 11 13 3 1 11 13 15
3 1 7 5 3 1 7 9
3 1 11 9 3 1 11 13
3 5 7 9 3 5 15 9
3 5 15 13 3 13 11 9
3 13 15 5 3 13 15 9
7 1 3 5 11 1 3 5
5 3 13 11 5 3 13 15
11 1 7 5 5 7 9 11
5 7 9 15 5 15 9 7
5 15 9 11 5 15 13 11
7 1 3 13 7 1 11 9
7 1 11 13 7 5 3 13
7 5 15 9 7 5 15 13
7 9 11 13 7 9 15 13
11 1 7 9 9 7 5 15
9 11 13 15 9 15 13 11
11 1 3 13 11 9 15 13
13 3 5 15 13 11 9 15
0 !NIMPHI
0 !NDON
0 !NACC
0 !NNB
0 !NGRP

View file

@ -5,3 +5,8 @@
2gly_EB-NEB.inp
#
2gly_CI-NEB.inp
# faster bonded correct gaussian and bug fix for stress tensor in shell model
2gly_B-NEB.inp
2gly_IT-NEB.inp
# faster bonded correct gaussian and bug fix for stress tensor in shell model
2gly_CI-NEB.inp

View file

@ -16,3 +16,9 @@
2gly_DIIS-OEP-NEB.inp
#
2gly_DIIS-SD-2.inp
# faster bonded correct gaussian and bug fix for stress tensor in shell model
2gly_DIIS-DNEB.inp
2gly_DIIS-SD-2.inp
2gly_DIIS-SD-NEB.inp
# faster bonded correct gaussian and bug fix for stress tensor in shell model
2gly_DIIS-NEB.inp

View file

@ -1 +1,3 @@
#
# faster bonded correct gaussian and bug fix for stress tensor in shell model
2gly_IT-NEB-res.inp