diff --git a/src/fist_force.F b/src/fist_force.F index e3bef9dcb4..9d192d32cf 100644 --- a/src/fist_force.F +++ b/src/fist_force.F @@ -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, & diff --git a/src/fist_nonbond_force.F b/src/fist_nonbond_force.F index a86ed43943..7f5d7d3492 100644 --- a/src/fist_nonbond_force.F +++ b/src/fist_nonbond_force.F @@ -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 diff --git a/tests/Fist/regtest-5/TEST_FILES_RESET b/tests/Fist/regtest-5/TEST_FILES_RESET index 4935f492bf..44b64c451a 100644 --- a/tests/Fist/regtest-5/TEST_FILES_RESET +++ b/tests/Fist/regtest-5/TEST_FILES_RESET @@ -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 diff --git a/tests/Fist/regtest-6/TEST_FILES b/tests/Fist/regtest-6/TEST_FILES index ac3d5e5111..1713bc11c8 100644 --- a/tests/Fist/regtest-6/TEST_FILES +++ b/tests/Fist/regtest-6/TEST_FILES @@ -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 diff --git a/tests/Fist/regtest-6/TEST_FILES_RESET b/tests/Fist/regtest-6/TEST_FILES_RESET index 4b016e2c7a..c9b18d7e54 100644 --- a/tests/Fist/regtest-6/TEST_FILES_RESET +++ b/tests/Fist/regtest-6/TEST_FILES_RESET @@ -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 diff --git a/tests/Fist/regtest-6/cubane_0.inp b/tests/Fist/regtest-6/cubane_0.inp new file mode 100644 index 0000000000..38b6ec2c8a --- /dev/null +++ b/tests/Fist/regtest-6/cubane_0.inp @@ -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 + + diff --git a/tests/Fist/regtest-6/cubane_15.inp b/tests/Fist/regtest-6/cubane_15.inp new file mode 100644 index 0000000000..3d3a0f8c2d --- /dev/null +++ b/tests/Fist/regtest-6/cubane_15.inp @@ -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 + + diff --git a/tests/Fist/regtest/TEST_FILES_RESET b/tests/Fist/regtest/TEST_FILES_RESET index f4d893e7ab..c3bc6ce3fd 100644 --- a/tests/Fist/regtest/TEST_FILES_RESET +++ b/tests/Fist/regtest/TEST_FILES_RESET @@ -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 diff --git a/tests/Fist/sample_pdb/cubane_0.pdb b/tests/Fist/sample_pdb/cubane_0.pdb new file mode 100644 index 0000000000..d70008cce6 --- /dev/null +++ b/tests/Fist/sample_pdb/cubane_0.pdb @@ -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 diff --git a/tests/Fist/sample_pdb/cubane_15.pdb b/tests/Fist/sample_pdb/cubane_15.pdb new file mode 100644 index 0000000000..623f1f4cf2 --- /dev/null +++ b/tests/Fist/sample_pdb/cubane_15.pdb @@ -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 diff --git a/tests/Fist/sample_pot/cubane.pot b/tests/Fist/sample_pot/cubane.pot new file mode 100644 index 0000000000..f5742c0dda --- /dev/null +++ b/tests/Fist/sample_pot/cubane.pot @@ -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 diff --git a/tests/Fist/sample_psf/cubane.psf b/tests/Fist/sample_psf/cubane.psf new file mode 100644 index 0000000000..d5313d7f49 --- /dev/null +++ b/tests/Fist/sample_psf/cubane.psf @@ -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 + diff --git a/tests/NEB/regtest-1/TEST_FILES_RESET b/tests/NEB/regtest-1/TEST_FILES_RESET index a46bd75534..0d989b62c2 100644 --- a/tests/NEB/regtest-1/TEST_FILES_RESET +++ b/tests/NEB/regtest-1/TEST_FILES_RESET @@ -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 diff --git a/tests/NEB/regtest-2/TEST_FILES_RESET b/tests/NEB/regtest-2/TEST_FILES_RESET index 1ed8aecd79..bbca0231ab 100644 --- a/tests/NEB/regtest-2/TEST_FILES_RESET +++ b/tests/NEB/regtest-2/TEST_FILES_RESET @@ -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 diff --git a/tests/NEB/regtest-3/TEST_FILES_RESET b/tests/NEB/regtest-3/TEST_FILES_RESET index 792d600548..814472ed3c 100644 --- a/tests/NEB/regtest-3/TEST_FILES_RESET +++ b/tests/NEB/regtest-3/TEST_FILES_RESET @@ -1 +1,3 @@ # +# faster bonded correct gaussian and bug fix for stress tensor in shell model +2gly_IT-NEB-res.inp