From 6aac6bb711b74b7e4585ad48baeabca4ec91fb3d Mon Sep 17 00:00:00 2001 From: Teodoro Laino Date: Fri, 7 Nov 2008 18:09:18 +0000 Subject: [PATCH] correct the analytical form of the dihedral potential svn-origin-rev: 7903 --- src/fist_intra_force.F | 42 ++++++----- src/force_field_types.F | 40 +++++----- src/force_fields.F | 2 +- src/force_fields_all.F | 28 +++---- src/force_fields_ext.F | 22 +++--- src/force_fields_input.F | 14 ++-- src/force_fields_util.F | 6 +- src/input_cp2k_mm.F | 8 +- src/mc_misc.F | 2 +- src/mol_force.F | 95 ++++++++---------------- src/topology_amber.F | 12 ++- tests/FE/regtest-1/TEST_FILES_RESET | 4 + tests/FE/regtest-2/TEST_FILES_RESET | 2 + tests/FE/regtest-3/TEST_FILES_RESET | 3 + tests/Fist/regtest-3/TEST_FILES_RESET | 4 + tests/Fist/regtest-4/TEST_FILES_RESET | 5 ++ tests/Fist/regtest-5/TEST_FILES_RESET | 9 +++ tests/Fist/regtest/TEST_FILES_RESET | 4 + tests/NEB/regtest-1/TEST_FILES_RESET | 5 ++ tests/NEB/regtest-2/TEST_FILES_RESET | 7 ++ tests/NEB/regtest-3/TEST_FILES_RESET | 3 + tests/NEB/regtest-4/TEST_FILES_RESET | 3 + tests/QMMM/QS/regtest-3/TEST_FILES_RESET | 3 + tests/QMMM/SE/regtest/TEST_FILES_RESET | 6 ++ 24 files changed, 178 insertions(+), 151 deletions(-) diff --git a/src/fist_intra_force.F b/src/fist_intra_force.F index 0ef508e648..e6589d47c6 100644 --- a/src/fist_intra_force.F +++ b/src/fist_intra_force.F @@ -85,12 +85,11 @@ CONTAINS index_a, index_b, index_c, index_d, ionef, ishell, itorsion, nbends, & nbonds, nimptors, nkind, nmol_per_kind, nonefs, nshell, ntorsions, nub REAL(KIND=dp) :: d12, d32, dist, dist1, dist2, energy, fscalar, id12, & - id32, is32, ism, isn, iw1, iw2, rab2, rinv, rinv2, s32, sm, sn, theta, & - w1, w2 + id32, is32, ism, isn, rab2, rinv, rinv2, s32, sm, sn, theta REAL(KIND=dp), DIMENSION(3) :: b12, b32, g1, g2, g3, gt1, & gt2, gt3, gt4, k1, k2, k3, & k4, rij, t12, t32, t34, t43, & - tm, tn, u1, u2 + tm, tn TYPE(bend_type), POINTER :: bend_list( : ) TYPE(bond_type), POINTER :: bond_list( : ) TYPE(cp_logger_type), POINTER :: logger @@ -310,26 +309,33 @@ CONTAINS index_d = torsion_list ( itorsion ) % d + first_atom - 1 t12 = particle_set ( index_a ) % r - particle_set ( index_b ) % r t32 = particle_set ( index_c ) % r - particle_set ( index_b ) % r + t34 = particle_set ( index_c ) % r - particle_set ( index_d ) % r t43 = particle_set ( index_d ) % r - particle_set ( index_c ) % r t12 = pbc(t12,cell) t32 = pbc(t32,cell) + t34 = pbc(t34,cell) t43 = pbc(t43,cell) - s32 = SQRT(DOT_PRODUCT(t32,t32)) + ! t12 x t32 + tm(1) = t12(2)*t32(3)-t32(2)*t12(3) + tm(2) = -t12(1)*t32(3)+t32(1)*t12(3) + tm(3) = t12(1)*t32(2)-t32(1)*t12(2) + ! t32 x t34 + tn(1) = t32(2)*t34(3)-t34(2)*t32(3) + tn(2) = -t32(1)*t34(3)+t34(1)*t32(3) + tn(3) = t32(1)*t34(2)-t34(1)*t32(2) + sm = SQRT ( DOT_PRODUCT(tm,tm) ) + ism = 1.0_dp/sm + sn = SQRT ( DOT_PRODUCT(tn,tn) ) + isn = 1.0_dp/sn + s32 = SQRT ( DOT_PRODUCT(t32,t32) ) is32 = 1.0_dp/s32 dist1 = DOT_PRODUCT(t12,t32) - dist2 = DOT_PRODUCT(t43,t32) - u1 = t12 - dist1*t32*is32**2 - u2 = t43 - dist2*t32*is32**2 - w1 = SQRT(DOT_PRODUCT(u1,u1)) - w2 = SQRT(DOT_PRODUCT(u2,u2)) - IF ((ABS(w1)0) THEN CALL reallocate(torsion_list(j)%torsion_kind%k,1,1) CALL reallocate(torsion_list(j)%torsion_kind%m,1,1) - CALL reallocate(torsion_list(j)%torsion_kind%cosphi0,1,1) + CALL reallocate(torsion_list(j)%torsion_kind%phi0,1,1) torsion_list(j)%torsion_kind%nmul = 1 torsion_list(j)%torsion_kind%m(1) = gro_info%torsion_m(itype) torsion_list(j)%torsion_kind%k(1) = gro_info%torsion_k(itype) - torsion_list(j)%torsion_kind%cosphi0(1) = gro_info%torsion_cosphi0(itype) + torsion_list(j)%torsion_kind%phi0(1) = gro_info%torsion_phi0(itype) ELSE CALL reallocate(torsion_list(j)%torsion_kind%k,1,1) CALL reallocate(torsion_list(j)%torsion_kind%m,1,1) - CALL reallocate(torsion_list(j)%torsion_kind%cosphi0,1,1) + CALL reallocate(torsion_list(j)%torsion_kind%phi0,1,1) torsion_list(j)%torsion_kind%nmul = 1 torsion_list(j)%torsion_kind%m(1) = gro_info%torsion_m(itype/k) torsion_list(j)%torsion_kind%k(1) = gro_info%torsion_k(itype/k) - torsion_list(j)%torsion_kind%cosphi0(1) = gro_info%torsion_cosphi0(itype/k) + torsion_list(j)%torsion_kind%phi0(1) = gro_info%torsion_phi0(itype/k) END IF torsion_list(j)%torsion_kind%id_type = gro_info%ff_gromos_type torsion_list(j)%id_type = gro_info%ff_gromos_type @@ -1343,11 +1343,11 @@ CONTAINS imul = torsion_list(j)%torsion_kind%nmul + 1 CALL reallocate(torsion_list(j)%torsion_kind%k,1,imul) CALL reallocate(torsion_list(j)%torsion_kind%m,1,imul) - CALL reallocate(torsion_list(j)%torsion_kind%cosphi0,1,imul) + CALL reallocate(torsion_list(j)%torsion_kind%phi0,1,imul) torsion_list(j)%torsion_kind%id_type = do_ff_charmm torsion_list(j)%torsion_kind%k(imul) = chm_info%torsion_k(k) torsion_list(j)%torsion_kind%m(imul) = chm_info%torsion_m(k) - torsion_list(j)%torsion_kind%cosphi0(imul) = chm_info%torsion_cosphi0(k) + torsion_list(j)%torsion_kind%phi0(imul) = chm_info%torsion_phi0(k) torsion_list(j)%torsion_kind%nmul = imul found = .TRUE. END IF @@ -1366,11 +1366,11 @@ CONTAINS imul = torsion_list(j)%torsion_kind%nmul + 1 CALL reallocate(torsion_list(j)%torsion_kind%k,1,imul) CALL reallocate(torsion_list(j)%torsion_kind%m,1,imul) - CALL reallocate(torsion_list(j)%torsion_kind%cosphi0,1,imul) + CALL reallocate(torsion_list(j)%torsion_kind%phi0,1,imul) torsion_list(j)%torsion_kind%id_type = do_ff_charmm torsion_list(j)%torsion_kind%k(imul) = chm_info%torsion_k(k) torsion_list(j)%torsion_kind%m(imul) = chm_info%torsion_m(k) - torsion_list(j)%torsion_kind%cosphi0(imul) = chm_info%torsion_cosphi0(k) + torsion_list(j)%torsion_kind%phi0(imul) = chm_info%torsion_phi0(k) torsion_list(j)%torsion_kind%nmul = imul found = .TRUE. END IF @@ -1392,11 +1392,11 @@ CONTAINS imul = torsion_list(j)%torsion_kind%nmul + 1 CALL reallocate(torsion_list(j)%torsion_kind%k,1,imul) CALL reallocate(torsion_list(j)%torsion_kind%m,1,imul) - CALL reallocate(torsion_list(j)%torsion_kind%cosphi0,1,imul) + CALL reallocate(torsion_list(j)%torsion_kind%phi0,1,imul) torsion_list(j)%torsion_kind%id_type = do_ff_amber torsion_list(j)%torsion_kind%k(imul) = amb_info%torsion_k(k) torsion_list(j)%torsion_kind%m(imul) = amb_info%torsion_m(k) - torsion_list(j)%torsion_kind%cosphi0(imul) = amb_info%torsion_cosphi0(k) + torsion_list(j)%torsion_kind%phi0(imul) = amb_info%torsion_phi0(k) torsion_list(j)%torsion_kind%nmul = imul found = .TRUE. END IF @@ -1415,11 +1415,11 @@ CONTAINS imul = torsion_list(j)%torsion_kind%nmul + 1 CALL reallocate(torsion_list(j)%torsion_kind%k,1,imul) CALL reallocate(torsion_list(j)%torsion_kind%m,1,imul) - CALL reallocate(torsion_list(j)%torsion_kind%cosphi0,1,imul) + CALL reallocate(torsion_list(j)%torsion_kind%phi0,1,imul) torsion_list(j)%torsion_kind%id_type = do_ff_amber torsion_list(j)%torsion_kind%k(imul) = amb_info%torsion_k(k) torsion_list(j)%torsion_kind%m(imul) = amb_info%torsion_m(k) - torsion_list(j)%torsion_kind%cosphi0(imul) = amb_info%torsion_cosphi0(k) + torsion_list(j)%torsion_kind%phi0(imul) = amb_info%torsion_phi0(k) torsion_list(j)%torsion_kind%nmul = imul found = .TRUE. END IF @@ -1441,11 +1441,11 @@ CONTAINS imul = torsion_list(j)%torsion_kind%nmul + 1 CALL reallocate(torsion_list(j)%torsion_kind%k,1,imul) CALL reallocate(torsion_list(j)%torsion_kind%m,1,imul) - CALL reallocate(torsion_list(j)%torsion_kind%cosphi0,1,imul) + CALL reallocate(torsion_list(j)%torsion_kind%phi0,1,imul) torsion_list(j)%torsion_kind%id_type = inp_info%torsion_kind(k) torsion_list(j)%torsion_kind%k(imul) = inp_info%torsion_k(k) torsion_list(j)%torsion_kind%m(imul) = inp_info%torsion_m(k) - torsion_list(j)%torsion_kind%cosphi0(imul) = COS(inp_info%torsion_cosphi0(k)) + torsion_list(j)%torsion_kind%phi0(imul) = inp_info%torsion_phi0(k) torsion_list(j)%torsion_kind%nmul = imul found = .TRUE. END IF diff --git a/src/force_fields_ext.F b/src/force_fields_ext.F index 7d11b7ab6b..57e9d0e2ee 100644 --- a/src/force_fields_ext.F +++ b/src/force_fields_ext.F @@ -81,7 +81,7 @@ CONTAINS itype, iw, jatom, ncon, & ntype, offset, stat LOGICAL :: failure, found - REAL(KIND=dp) :: cost2, csq, sdet + REAL(KIND=dp) :: cosphi0, cost2, csq, sdet TYPE(cp_logger_type), POINTER :: logger TYPE(cp_parser_type), POINTER :: parser TYPE(gromos_info_type), POINTER :: gro_info @@ -236,12 +236,13 @@ CONTAINS CALL parser_get_object(parser,ntype,error=error) CALL reallocate(gro_info%torsion_k,1,ntype) CALL reallocate(gro_info%torsion_m,1,ntype) - CALL reallocate(gro_info%torsion_cosphi0,1,ntype) + CALL reallocate(gro_info%torsion_phi0,1,ntype) DO itype=1,ntype CALL parser_get_next_line(parser,1,error=error) CALL parser_get_object(parser,gro_info%torsion_k(itype),error=error) - CALL parser_get_object(parser,gro_info%torsion_cosphi0(itype),error=error) + CALL parser_get_object(parser,cosphi0,error=error) CALL parser_get_object(parser,gro_info%torsion_m(itype),error=error) + gro_info%torsion_phi0(itype) = ACOS(cosphi0) gro_info%torsion_k(itype) = cp_unit_to_cp2k(gro_info%torsion_k(itype),"kjmol",error=error) IF(iw>0) WRITE(iw,*) "GROMOS_FF| PUT DIHEDRALTYPE INFO HERE!!!!" END DO @@ -523,14 +524,14 @@ CONTAINS CALL reallocate(chm_info%torsion_d,1,ntorsion) CALL reallocate(chm_info%torsion_k,1,ntorsion) CALL reallocate(chm_info%torsion_m,1,ntorsion) - CALL reallocate(chm_info%torsion_cosphi0,1,ntorsion) + CALL reallocate(chm_info%torsion_phi0,1,ntorsion) chm_info%torsion_a(ntorsion) = string chm_info%torsion_b(ntorsion) = string2 chm_info%torsion_c(ntorsion) = string3 chm_info%torsion_d(ntorsion) = string4 CALL parser_get_object(parser,chm_info%torsion_k(ntorsion),error=error) CALL parser_get_object(parser,chm_info%torsion_m(ntorsion),error=error) - CALL parser_get_object(parser,chm_info%torsion_cosphi0(ntorsion),error=error) + CALL parser_get_object(parser,chm_info%torsion_phi0(ntorsion),error=error) IF(iw>0) WRITE(iw,*) " CHM TORSION ",ntorsion,& TRIM(chm_info%torsion_a(ntorsion))," ",& TRIM(chm_info%torsion_b(ntorsion))," ",& @@ -538,14 +539,12 @@ CONTAINS TRIM(chm_info%torsion_d(ntorsion))," ",& chm_info%torsion_k(ntorsion),& chm_info%torsion_m(ntorsion),& - chm_info%torsion_cosphi0(ntorsion) + chm_info%torsion_phi0(ntorsion) ! Do some units conversion into internal atomic units - chm_info%torsion_cosphi0(ntorsion) = cp_unit_to_cp2k(chm_info%torsion_cosphi0(ntorsion),& + chm_info%torsion_phi0(ntorsion) = cp_unit_to_cp2k(chm_info%torsion_phi0(ntorsion),& "deg",error=error) chm_info%torsion_k(ntorsion) = cp_unit_to_cp2k(chm_info%torsion_k(ntorsion),"kcalmol",& error=error) - chm_info%torsion_cosphi0(ntorsion)=& - COS(chm_info%torsion_cosphi0(ntorsion)) CALL charmm_get_next_line(parser,1,error=error) END DO ELSE @@ -764,11 +763,10 @@ CONTAINS TRIM(amb_info%torsion_d(i))," ",& amb_info%torsion_k(i),& amb_info%torsion_m(i),& - amb_info%torsion_cosphi0(i) + amb_info%torsion_phi0(i) ! Do some units conversion into internal atomic units - ! For AMBER forcefield the cosphi0 contains the phi0 since the potential - amb_info%torsion_cosphi0(i) = cp_unit_to_cp2k(amb_info%torsion_cosphi0(i),"rad",error=error) + amb_info%torsion_phi0(i) = cp_unit_to_cp2k(amb_info%torsion_phi0(i),"rad",error=error) amb_info%torsion_k(i) = cp_unit_to_cp2k(amb_info%torsion_k(i),"kcalmol",error=error) END DO diff --git a/src/force_fields_input.F b/src/force_fields_input.F index 1fb444af82..86858260eb 100644 --- a/src/force_fields_input.F +++ b/src/force_fields_input.F @@ -336,10 +336,10 @@ CPSourceFileRef,& CALL reallocate(inp_info%torsion_d,1,ntors) CALL reallocate(inp_info%torsion_k,1,ntors) CALL reallocate(inp_info%torsion_m,1,ntors) - CALL reallocate(inp_info%torsion_cosphi0,1,ntors) - CALL read_torsions_section(inp_info%torsion_kind,inp_info%torsion_a, inp_info%torsion_b, inp_info%torsion_c,& - inp_info%torsion_d,inp_info%torsion_k,inp_info%torsion_cosphi0, inp_info%torsion_m,& - tmp_section, ntot, error ) + CALL reallocate(inp_info%torsion_phi0,1,ntors) + CALL read_torsions_section(inp_info%torsion_kind,inp_info%torsion_a, inp_info%torsion_b,& + inp_info%torsion_c,inp_info%torsion_d,inp_info%torsion_k,inp_info%torsion_phi0,& + inp_info%torsion_m,tmp_section, ntot, error ) END IF END IF @@ -1209,12 +1209,12 @@ CPSourceFileRef,& !> \author teo ! ***************************************************************************** SUBROUTINE read_torsions_section(torsion_kind,torsion_a, torsion_b, torsion_c, torsion_d, torsion_k,& - torsion_cosphi0, torsion_m, section, start, error ) + torsion_phi0, torsion_m, section, start, error ) INTEGER, DIMENSION(:), POINTER :: torsion_kind CHARACTER(LEN=default_string_length), & DIMENSION(:), POINTER :: torsion_a, torsion_b, & torsion_c, torsion_d - REAL(KIND=dp), DIMENSION(:), POINTER :: torsion_k, torsion_cosphi0 + REAL(KIND=dp), DIMENSION(:), POINTER :: torsion_k, torsion_phi0 INTEGER, DIMENSION(:), POINTER :: torsion_m TYPE(section_vals_type), POINTER :: section INTEGER, INTENT(IN) :: start @@ -1243,7 +1243,7 @@ CPSourceFileRef,& CALL uppercase(torsion_c(start+isec)) CALL uppercase(torsion_d(start+isec)) CALL section_vals_val_get(section,"K",i_rep_section=isec,r_val=torsion_k(start+isec),error=error) - CALL section_vals_val_get(section,"PHI0",i_rep_section=isec,r_val=torsion_cosphi0(start+isec),error=error) + CALL section_vals_val_get(section,"PHI0",i_rep_section=isec,r_val=torsion_phi0(start+isec),error=error) CALL section_vals_val_get(section,"M",i_rep_section=isec,i_val=torsion_m(start+isec),error=error) END DO END SUBROUTINE read_torsions_section diff --git a/src/force_fields_util.F b/src/force_fields_util.F index c6191a3f91..0c4304b6e1 100644 --- a/src/force_fields_util.F +++ b/src/force_fields_util.F @@ -976,16 +976,16 @@ CONTAINS new_torsion_kind_set(counter)= torsion_kind_set(itorsion) i = SIZE(torsion_kind_set(itorsion)%m) j = SIZE(torsion_kind_set(itorsion)%k) - k = SIZE(torsion_kind_set(itorsion)%cosphi0) + k = SIZE(torsion_kind_set(itorsion)%phi0) ALLOCATE(new_torsion_kind_set(counter)%m(i),STAT=stat) CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) ALLOCATE(new_torsion_kind_set(counter)%k(i),STAT=stat) CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) - ALLOCATE(new_torsion_kind_set(counter)%cosphi0(i),STAT=stat) + ALLOCATE(new_torsion_kind_set(counter)%phi0(i),STAT=stat) CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) new_torsion_kind_set(counter)%m= torsion_kind_set(itorsion)%m new_torsion_kind_set(counter)%k= torsion_kind_set(itorsion)%k - new_torsion_kind_set(counter)%cosphi0= torsion_kind_set(itorsion)%cosphi0 + new_torsion_kind_set(counter)%phi0= torsion_kind_set(itorsion)%phi0 END IF END DO counter=0 diff --git a/src/input_cp2k_mm.F b/src/input_cp2k_mm.F index 1f3c576397..9e1e918537 100644 --- a/src/input_cp2k_mm.F +++ b/src/input_cp2k_mm.F @@ -490,10 +490,10 @@ CONTAINS description="Define the kind of torsion potential",& usage="KIND CHARMM",& enum_c_vals=s2a("CHARMM","G87","G96","AMBER"),& - enum_desc=s2a("Functional Form (CHARMM|G87|G96): K * [ 1 + cos[M*PHI] * cos[PHI0]]",& - "Functional Form (CHARMM|G87|G96): K * [ 1 + cos[M*PHI] * cos[PHI0]]",& - "Functional Form (CHARMM|G87|G96): K * [ 1 + cos[M*PHI] * cos[PHI0]]",& - "Function Form AMBER : K * [ 1 + cos[M*PHI - PHI0]]"),& + enum_desc=s2a("Functional Form (CHARMM|G87|G96|AMBER): K * [ 1 + cos[M*PHI - PHI0]]",& + "Functional Form (CHARMM|G87|G96|AMBER): K * [ 1 + cos[M*PHI - PHI0]]",& + "Functional Form (CHARMM|G87|G96|AMBER): K * [ 1 + cos[M*PHI - PHI0]]",& + "Functional Form (CHARMM|G87|G96|AMBER): K * [ 1 + cos[M*PHI - PHI0]]"),& enum_i_vals=(/ do_ff_charmm,& do_ff_g87,& do_ff_g96,& diff --git a/src/mc_misc.F b/src/mc_misc.F index 578431179b..4ca1d66ee2 100644 --- a/src/mc_misc.F +++ b/src/mc_misc.F @@ -767,7 +767,7 @@ SUBROUTINE mc_make_dat_file(coordinates,natoms_tot,box_length,filename,& TRIM(ADJUSTL(name_c)),TRIM(ADJUSTL(name_d)) WRITE(unit,20003) " M ",torsion_list(itorsion)%torsion_kind%m(imul) WRITE(unit,20002) " K ",torsion_list(itorsion)%torsion_kind%k(imul) - WRITE(unit,20002) " PHI0 ",torsion_list(itorsion)%torsion_kind%cosphi0(imul) + WRITE(unit,20002) " PHI0 ",torsion_list(itorsion)%torsion_kind%phi0(imul) WRITE(unit,10000) " &END TORSION" ENDDO diff --git a/src/mol_force.F b/src/mol_force.F index 5a497f2061..180beebb72 100644 --- a/src/mol_force.F +++ b/src/mol_force.F @@ -191,12 +191,13 @@ CONTAINS !> none !> \author DG ! ***************************************************************************** - RECURSIVE SUBROUTINE force_torsions (id_type, is32, dist1, dist2, u1, u2, iw1, iw2, & - k, cosphi0, m, gt1, gt2, gt3, gt4, energy, fscalar ) + SUBROUTINE force_torsions (id_type, s32, is32, ism, isn, dist1, dist2, tm,& + tn, t12, t32, k, phi0, m, gt1, gt2, gt3, gt4, energy, fscalar ) INTEGER, INTENT(IN) :: id_type - REAL(KIND=dp), INTENT(IN) :: is32, dist1, dist2 - REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: u1, u2 - REAL(KIND=dp), INTENT(IN) :: iw1, iw2, k, cosphi0 + REAL(KIND=dp), INTENT(IN) :: s32, is32, ism, isn, dist1, & + dist2 + REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: tm, tn, t12, t32 + REAL(KIND=dp), INTENT(IN) :: k, phi0 INTEGER, INTENT(IN) :: m REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: gt1, gt2, gt3, gt4 REAL(KIND=dp), INTENT(OUT) :: energy, fscalar @@ -204,66 +205,32 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'force_torsions', & routineP = moduleN//':'//routineN - REAL(KIND=dp) :: cosphi, lcosphi0, trig + LOGICAL :: failure + REAL(KIND=dp) :: cosphi, phi - cosphi = DOT_PRODUCT ( u1, u2 ) * iw1 * iw2 + failure = .FALSE. + cosphi = DOT_PRODUCT ( tm, tn ) * ism * isn + IF(cosphi > 1.0_dp) cosphi = 1.0_dp + IF(cosphi < -1.0_dp) cosphi = -1.0_dp + phi = SIGN( ACOS ( cosphi ), DOT_PRODUCT ( t12, tn )) - ! compute energy + ! Select force field SELECT CASE (id_type) - CASE (do_ff_charmm,do_ff_g87,do_ff_g96) - SELECT CASE ( m ) - CASE ( 0 ) - trig = 1.0_dp - fscalar = 0.0_dp - CASE ( 1 ) - trig = cosphi - fscalar = 1.0_dp - CASE ( 2 ) - trig = 2.0_dp * cosphi* cosphi - 1.0_dp - fscalar = 4.0_dp * cosphi - CASE ( 3 ) - trig = 4.0_dp * cosphi * cosphi * cosphi - & - 3.0_dp * cosphi - fscalar = 12.0_dp * cosphi * cosphi - 3.0_dp - CASE ( 4 ) - trig = 8.0_dp * cosphi * cosphi * cosphi * cosphi - & - 8.0_dp * cosphi * cosphi + 1.0_dp - fscalar = 32.0_dp * cosphi * cosphi * cosphi - 16.0_dp * cosphi - CASE ( 5 ) - trig = 16.0_dp * cosphi * cosphi * cosphi * cosphi * cosphi - & - 20.0_dp * cosphi * cosphi * cosphi + & - 5.0_dp * cosphi - fscalar = 80.0_dp * cosphi * cosphi * cosphi * cosphi - & - 60.0_dp * cosphi * cosphi + 5.0_dp - CASE ( 6 ) - trig = 32.0_dp * cosphi * cosphi * cosphi * cosphi * cosphi * cosphi - & - 48.0_dp * cosphi * cosphi * cosphi * cosphi + & - 18.0_dp * cosphi * cosphi - 1.0_dp - fscalar = 192.0_dp * cosphi * cosphi * cosphi * cosphi * cosphi - & - 192.0_dp * cosphi * cosphi * cosphi + 36.0_dp * cosphi - CASE DEFAULT - CALL stop_program (routineP,"torsion parameter too large") - END SELECT - energy = k * ( 1.0_dp + cosphi0 * trig ) + CASE (do_ff_charmm,do_ff_g87,do_ff_g96, do_ff_amber) + ! compute energy + energy = k * ( 1.0_dp + COS(m*phi-phi0)) ! compute fscalar - fscalar = -fscalar * cosphi0 * k - - ! compute the gradients - gt1 = ( u2 * iw2 - cosphi * u1 * iw1 ) * iw1 - gt4 = ( u1 * iw1 - cosphi * u2 * iw2 ) * iw2 - gt2 = ( dist1 * is32 ** 2 - 1.0_dp ) * gt1 + dist2 * is32 ** 2 * gt4 - gt3 = ( - dist2 * is32 ** 2 - 1.0_dp ) * gt4 - dist1 * is32 ** 2 * gt1 - - CASE (do_ff_amber) - ! Use the same expression of CHARMM for AMBER (to be checked!) - lcosphi0 = COS(cosphi0) - CALL force_torsions (do_ff_charmm, is32, dist1, dist2, u1, u2, iw1, iw2, & - k, lcosphi0, m, gt1, gt2, gt3, gt4, energy, fscalar ) + fscalar = k * m * SIN(m*phi-phi0) CASE DEFAULT CALL stop_program (routineP,"unmatched torsion kind") END SELECT + ! compute the gradients + gt1 = ( s32 * ism * ism ) * tm + gt4 = -( s32 * isn * isn ) * tn + gt2 = ( dist1 * is32 ** 2 - 1.0_dp ) * gt1 - dist2 * is32 ** 2 * gt4 + gt3 = ( dist2 * is32 ** 2 - 1.0_dp ) * gt4 - dist1 * is32 ** 2 * gt1 END SUBROUTINE force_torsions ! ***************************************************************************** @@ -286,14 +253,15 @@ CONTAINS routineP = moduleN//':'//routineN REAL(KIND=dp), PARAMETER :: f12 = 1.0_dp/2.0_dp + LOGICAL :: failure REAL(KIND=dp) :: cosphi, phi -! define cosphi - + failure = .FALSE. + ! define cosphi cosphi = DOT_PRODUCT ( tm, tn ) * ism * isn - IF(cosphi > 1.0_dp) cosphi = 1.0_dp + IF(cosphi > 1.0_dp) cosphi = 1.0_dp IF(cosphi < -1.0_dp) cosphi = -1.0_dp - phi = SIGN ( ACOS ( cosphi ), DOT_PRODUCT ( t12, tn ) ) + phi = SIGN( ACOS ( cosphi ), DOT_PRODUCT ( t12, tn )) SELECT CASE (id_type) CASE (do_ff_charmm) @@ -315,11 +283,10 @@ CONTAINS END SELECT ! compute the gradients - gt1 = ( s32* ism * ism ) * tm + gt1 = ( s32 * ism * ism ) * tm gt4 = -( s32 * isn * isn ) * tn - gt2 = ( dist1 * is32 ** 2 - 1.0_dp ) * gt1 - dist2 * is32 ** 2 * gt4 - gt3 = ( dist2 * is32 ** 2 - 1.0_dp ) * gt4 - dist1 * is32 ** 2 * gt1 - + gt2 = ( dist1 * is32 ** 2 - 1.0_dp ) * gt1 - dist2 * is32 ** 2 * gt4 + gt3 = ( dist2 * is32 ** 2 - 1.0_dp ) * gt4 - dist1 * is32 ** 2 * gt1 END SUBROUTINE force_imp_torsions ! ***************************************************************************** diff --git a/src/topology_amber.F b/src/topology_amber.F index 1adc5b5adb..9cee74196e 100644 --- a/src/topology_amber.F +++ b/src/topology_amber.F @@ -854,9 +854,7 @@ CPSourceFileRef,& ! ---------------------------------------------------------- ! Force Fields informations related to torsions - ! in amb_info%cosphi0 we store PHI0 and not its cosine, this - ! is just to avoid an additional field in the torsional kind. - ! Just beware of it! + ! in amb_info%phi0 we store PHI0 ! ---------------------------------------------------------- CALL reallocate(amb_info%torsion_a, 1,buffer_size) CALL reallocate(amb_info%torsion_b, 1,buffer_size) @@ -864,17 +862,17 @@ CPSourceFileRef,& CALL reallocate(amb_info%torsion_d, 1,buffer_size) CALL reallocate(amb_info%torsion_k, 1,buffer_size) CALL reallocate(amb_info%torsion_m, 1,buffer_size) - CALL reallocate(amb_info%torsion_cosphi0, 1,buffer_size) + CALL reallocate(amb_info%torsion_phi0, 1,buffer_size) nsize = 0 ! Torsions containing hydrogens CALL post_process_torsions_info(amb_info%torsion_a, amb_info%torsion_b,& amb_info%torsion_c, amb_info%torsion_d, amb_info%torsion_k, & - amb_info%torsion_m, amb_info%torsion_cosphi0, particle_set, nsize,& + amb_info%torsion_m, amb_info%torsion_phi0, particle_set, nsize,& nphih, iph, jph, kph, lph, icph, pk, pn, phase, error) ! Torsions non-containing hydrogens CALL post_process_torsions_info(amb_info%torsion_a, amb_info%torsion_b,& amb_info%torsion_c, amb_info%torsion_d, amb_info%torsion_k, & - amb_info%torsion_m, amb_info%torsion_cosphi0, particle_set, nsize,& + amb_info%torsion_m, amb_info%torsion_phi0, particle_set, nsize,& nphia, ip, jp, kp, lp, icp, pk, pn, phase, error) ! Shrink arrays size to the minimal request CALL reallocate(amb_info%torsion_a, 1,nsize) @@ -883,7 +881,7 @@ CPSourceFileRef,& CALL reallocate(amb_info%torsion_d, 1,nsize) CALL reallocate(amb_info%torsion_k, 1,nsize) CALL reallocate(amb_info%torsion_m, 1,nsize) - CALL reallocate(amb_info%torsion_cosphi0, 1,nsize) + CALL reallocate(amb_info%torsion_phi0, 1,nsize) ! ---------------------------------------------------------- ! Post process of LJ parameters diff --git a/tests/FE/regtest-1/TEST_FILES_RESET b/tests/FE/regtest-1/TEST_FILES_RESET index 6c30bccd8b..dd2f393bd6 100644 --- a/tests/FE/regtest-1/TEST_FILES_RESET +++ b/tests/FE/regtest-1/TEST_FILES_RESET @@ -14,3 +14,7 @@ Solv_alch_chng_simpl.inp Solv_alch_chng.inp Solv_alch_chng_dist.inp Solv_alch_chng_simpl.inp +# numerics +Solv_alch_chng_dist.inp +Solv_alch_chng.inp +Solv_alch_chng_simpl.inp diff --git a/tests/FE/regtest-2/TEST_FILES_RESET b/tests/FE/regtest-2/TEST_FILES_RESET index fbf895ee15..06033e1ccf 100644 --- a/tests/FE/regtest-2/TEST_FILES_RESET +++ b/tests/FE/regtest-2/TEST_FILES_RESET @@ -6,3 +6,5 @@ Solv_alch_chng.inp Solv_alch_chng.inp # eps0 Solv_alch_chng.inp +# numerics +Solv_alch_chng.inp diff --git a/tests/FE/regtest-3/TEST_FILES_RESET b/tests/FE/regtest-3/TEST_FILES_RESET index a3a9d5b6f0..81a61a8bd7 100644 --- a/tests/FE/regtest-3/TEST_FILES_RESET +++ b/tests/FE/regtest-3/TEST_FILES_RESET @@ -10,3 +10,6 @@ Solv_alch_chng_res.inp # eps0 Solv_alch_chng.inp Solv_alch_chng_res.inp +# numerics +Solv_alch_chng.inp +Solv_alch_chng_res.inp diff --git a/tests/Fist/regtest-3/TEST_FILES_RESET b/tests/Fist/regtest-3/TEST_FILES_RESET index f09fd95176..401a17fdcd 100644 --- a/tests/Fist/regtest-3/TEST_FILES_RESET +++ b/tests/Fist/regtest-3/TEST_FILES_RESET @@ -99,3 +99,7 @@ nh3_dist.inp nh3_dist_npt.inp NaCl.inp NaCl-H2O.inp +# numerics +ace_ala_nme-ambconn.inp +ace_ala_nme-amber.inp +ace_ala_nme.inp diff --git a/tests/Fist/regtest-4/TEST_FILES_RESET b/tests/Fist/regtest-4/TEST_FILES_RESET index d1786de365..501b87f67b 100644 --- a/tests/Fist/regtest-4/TEST_FILES_RESET +++ b/tests/Fist/regtest-4/TEST_FILES_RESET @@ -164,3 +164,8 @@ nh3-meta-walks_2.inp nh3-meta-walks_2r.inp # yet numerics nh3-meta-walks_2r.inp +# numerics +cycbut.inp +cyhex.inp +gly_amber.inp +loop.inp diff --git a/tests/Fist/regtest-5/TEST_FILES_RESET b/tests/Fist/regtest-5/TEST_FILES_RESET index 016eaf557d..0aba82ea29 100644 --- a/tests/Fist/regtest-5/TEST_FILES_RESET +++ b/tests/Fist/regtest-5/TEST_FILES_RESET @@ -158,3 +158,12 @@ JAC_gen.inp # numerics after new units converter N3dye_vib_bfgs2.inp N3dye_vib_bfgs.inp +# numerics +arginine0.inp +arginine_crossC.inp +arginine.inp +N3dye_geoopt.inp +N3dye_vib_bfgs2.inp +N3dye_vib_bfgs.inp +N3dye_vib_restart_vec.inp +phenobenz.inp diff --git a/tests/Fist/regtest/TEST_FILES_RESET b/tests/Fist/regtest/TEST_FILES_RESET index f482ecf5af..d23b372167 100644 --- a/tests/Fist/regtest/TEST_FILES_RESET +++ b/tests/Fist/regtest/TEST_FILES_RESET @@ -463,3 +463,7 @@ multipole_quadrupole.dbg_f_rec.inp multipole_dipole.dbg_f_rec.inp multipole_dip_qu.dbg_f_rec.inp multipole_quadrupole.dbg_f_rec.inp +# numerics +multi_frag_check.inp +multi_frag.inp +multi_psf.inp diff --git a/tests/NEB/regtest-1/TEST_FILES_RESET b/tests/NEB/regtest-1/TEST_FILES_RESET index a183198df2..856d809d51 100644 --- a/tests/NEB/regtest-1/TEST_FILES_RESET +++ b/tests/NEB/regtest-1/TEST_FILES_RESET @@ -44,3 +44,8 @@ # numerics after new units converter 2gly_CI-NEB.inp 2gly_IT-NEB.inp +# numerics +2gly_B-NEB.inp +2gly_CI-NEB.inp +2gly_EB-NEB.inp +2gly_IT-NEB.inp diff --git a/tests/NEB/regtest-2/TEST_FILES_RESET b/tests/NEB/regtest-2/TEST_FILES_RESET index af0950127b..f85e3eb58d 100644 --- a/tests/NEB/regtest-2/TEST_FILES_RESET +++ b/tests/NEB/regtest-2/TEST_FILES_RESET @@ -58,3 +58,10 @@ 2gly_DIIS-SD-2.inp # numerics 2gly_DIIS-SD-NEB.inp +# numerics +2gly_DIIS-DNEB.inp +2gly_DIIS-NEB.inp +2gly_DIIS-OEP-NEB.inp +2gly_DIIS-SD-2.inp +2gly_DIIS-SD-NEB.inp +2gly_DIIS-SM.inp diff --git a/tests/NEB/regtest-3/TEST_FILES_RESET b/tests/NEB/regtest-3/TEST_FILES_RESET index f39cb4a938..a7ee9f2e5b 100644 --- a/tests/NEB/regtest-3/TEST_FILES_RESET +++ b/tests/NEB/regtest-3/TEST_FILES_RESET @@ -13,3 +13,6 @@ 2gly_IT-NEB-res.inp # numerics 2gly_IT-NEB-res.inp +# numerics +2gly_IT-NEB.inp +2gly_IT-NEB-res.inp diff --git a/tests/NEB/regtest-4/TEST_FILES_RESET b/tests/NEB/regtest-4/TEST_FILES_RESET index fde1352309..a783c8c652 100644 --- a/tests/NEB/regtest-4/TEST_FILES_RESET +++ b/tests/NEB/regtest-4/TEST_FILES_RESET @@ -5,3 +5,6 @@ # numerics after new units converter 2gly_IT-NEB-CV.inp 2gly_IT-NEB-CV-res.inp +# numerics +2gly_IT-NEB-CV.inp +2gly_IT-NEB-CV-res.inp diff --git a/tests/QMMM/QS/regtest-3/TEST_FILES_RESET b/tests/QMMM/QS/regtest-3/TEST_FILES_RESET index 390e86d8cf..320959a2e1 100644 --- a/tests/QMMM/QS/regtest-3/TEST_FILES_RESET +++ b/tests/QMMM/QS/regtest-3/TEST_FILES_RESET @@ -413,3 +413,6 @@ constr_hb_mm.inp constr_hb_qm.inp # corrected build of ppnl lists C4H10-qmmm-gauss-0.inp +# numerics +constr_hb_mm.inp +constr_hb_qm.inp diff --git a/tests/QMMM/SE/regtest/TEST_FILES_RESET b/tests/QMMM/SE/regtest/TEST_FILES_RESET index 1e7048b40d..35830ca183 100644 --- a/tests/QMMM/SE/regtest/TEST_FILES_RESET +++ b/tests/QMMM/SE/regtest/TEST_FILES_RESET @@ -96,3 +96,9 @@ water_3_fdbg_inv.inp water_3_walls.inp water_3_walls_q.inp zwitt.inp +# numerics +mol_CSVR_gen1.inp +mol_CSVR_gen2.inp +mol_CSVR_gen3.inp +mol.inp +mol_sph_cut.inp