diff --git a/src/force_field_types.F b/src/force_field_types.F index feac863ab5..d682ad055a 100644 --- a/src/force_field_types.F +++ b/src/force_field_types.F @@ -153,6 +153,10 @@ MODULE force_field_types CHARACTER(LEN=default_string_length), POINTER :: nonbond_a(:) REAL(KIND=dp), POINTER :: nonbond_eps(:) REAL(KIND=dp), POINTER :: nonbond_rmin2(:) + INTEGER, POINTER :: raw_torsion_id(:, :) + REAL(KIND=dp), POINTER :: raw_torsion_k(:) + REAL(KIND=dp), POINTER :: raw_torsion_m(:) + REAL(KIND=dp), POINTER :: raw_torsion_phi0(:) END TYPE amber_info_type ! ************************************************************************************************** @@ -360,6 +364,10 @@ CONTAINS NULLIFY (amb_info%nonbond_a) NULLIFY (amb_info%nonbond_eps) NULLIFY (amb_info%nonbond_rmin2) + NULLIFY (amb_info%raw_torsion_id) + NULLIFY (amb_info%raw_torsion_k) + NULLIFY (amb_info%raw_torsion_m) + NULLIFY (amb_info%raw_torsion_phi0) END SUBROUTINE init_amber_info @@ -720,6 +728,18 @@ CONTAINS IF (ASSOCIATED(amb_info%nonbond_rmin2)) THEN DEALLOCATE (amb_info%nonbond_rmin2) END IF + IF (ASSOCIATED(amb_info%raw_torsion_id)) THEN + DEALLOCATE (amb_info%raw_torsion_id) + END IF + IF (ASSOCIATED(amb_info%raw_torsion_k)) THEN + DEALLOCATE (amb_info%raw_torsion_k) + END IF + IF (ASSOCIATED(amb_info%raw_torsion_m)) THEN + DEALLOCATE (amb_info%raw_torsion_m) + END IF + IF (ASSOCIATED(amb_info%raw_torsion_phi0)) THEN + DEALLOCATE (amb_info%raw_torsion_phi0) + END IF END SUBROUTINE deallocate_amb_info diff --git a/src/force_fields_all.F b/src/force_fields_all.F index b5d104e41b..bb8e027c1d 100644 --- a/src/force_fields_all.F +++ b/src/force_fields_all.F @@ -439,7 +439,7 @@ CONTAINS ntorsion INTEGER, DIMENSION(:), POINTER :: molecule_list INTEGER, POINTER :: map_torsion_kind(:) - LOGICAL :: found + LOGICAL :: chk_reverse, found TYPE(atomic_kind_type), POINTER :: atomic_kind TYPE(molecule_kind_type), POINTER :: molecule_kind TYPE(molecule_type), POINTER :: molecule @@ -447,6 +447,11 @@ CONTAINS TYPE(torsion_type), DIMENSION(:), POINTER :: torsion_list CALL timeset(routineN, handle2) + + ! Now decide whether we need to check D-C-B-A type combination in addtion to usual A-B-C-D + ! We don't need it for Amber FF + chk_reverse = (ff_type%ff_type /= do_ff_amber) + DO i = 1, SIZE(molecule_kind_set) molecule_kind => molecule_kind_set(i) CALL get_molecule_kind(molecule_kind=molecule_kind, & @@ -503,7 +508,8 @@ CONTAINS ((name_atm_b) == (name_atm_b2)) .AND. & ((name_atm_c) == (name_atm_c2)) .AND. & ((name_atm_d) == (name_atm_d2))) .OR. & - (((name_atm_a) == (name_atm_d2)) .AND. & + (chk_reverse .AND. & + ((name_atm_a) == (name_atm_d2)) .AND. & ((name_atm_b) == (name_atm_c2)) .AND. & ((name_atm_c) == (name_atm_b2)) .AND. & ((name_atm_d) == (name_atm_a2)))) THEN @@ -1263,8 +1269,10 @@ CONTAINS CHARACTER(LEN=default_string_length) :: ldum, name_atm_a, name_atm_b, & name_atm_c, name_atm_d INTEGER :: atm_a, atm_b, atm_c, atm_d, first, & - handle2, i, imul, itype, j, k, last, & - natom, ntorsion + handle2, i, imul, itype, j, k, k_end, & + k_start, last, natom, ntorsion, & + raw_parm_id + INTEGER, DIMENSION(4) :: glob_atm_id INTEGER, DIMENSION(:), POINTER :: molecule_list LOGICAL :: found, only_qm TYPE(atomic_kind_type), POINTER :: atomic_kind @@ -1384,52 +1392,43 @@ CONTAINS END IF ! loop over params from AMBER + ! Assign real parameters from Amber PRMTOP file using global atom indices + ! Type-based assignment is prone to errors IF (ASSOCIATED(amb_info%torsion_a)) THEN - DO k = 1, SIZE(amb_info%torsion_a) - IF ((((amb_info%torsion_a(k)) == (name_atm_a)) .AND. & - ((amb_info%torsion_b(k)) == (name_atm_b)) .AND. & - ((amb_info%torsion_c(k)) == (name_atm_c)) .AND. & - ((amb_info%torsion_d(k)) == (name_atm_d))) .OR. & - (((amb_info%torsion_a(k)) == (name_atm_d)) .AND. & - ((amb_info%torsion_b(k)) == (name_atm_c)) .AND. & - ((amb_info%torsion_c(k)) == (name_atm_b)) .AND. & - ((amb_info%torsion_d(k)) == (name_atm_a)))) THEN + ! Get global atom indices + glob_atm_id(1) = atm_a + first - 1 + glob_atm_id(2) = atm_b + first - 1 + glob_atm_id(3) = atm_c + first - 1 + glob_atm_id(4) = atm_d + first - 1 + + ! Search sorted array of raw torsion parameters + ! The array can be too long for linear lookup + ! Use binary search for first atom index + k_start = bsearch_leftmost_2d(amb_info%raw_torsion_id, glob_atm_id(1)) + k_end = UBOUND(amb_info%raw_torsion_id, DIM=2) + + ! If not found, skip the loop + IF (k_start /= 0) THEN + + DO k = k_start, k_end + IF (glob_atm_id(1) < amb_info%raw_torsion_id(1, k)) EXIT + IF (ANY((glob_atm_id - amb_info%raw_torsion_id(1:4, k)) /= 0)) CYCLE + + raw_parm_id = amb_info%raw_torsion_id(5, k) 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%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%phi0(imul) = amb_info%torsion_phi0(k) + torsion_list(j)%torsion_kind%k(imul) = amb_info%raw_torsion_k(raw_parm_id) + torsion_list(j)%torsion_kind%m(imul) = NINT(amb_info%raw_torsion_m(raw_parm_id)) + torsion_list(j)%torsion_kind%phi0(imul) = amb_info%raw_torsion_phi0(raw_parm_id) torsion_list(j)%torsion_kind%nmul = imul found = .TRUE. - END IF - END DO - - IF (.NOT. found) THEN - DO k = 1, SIZE(amb_info%torsion_a) - IF ((((amb_info%torsion_a(k)) == ("X")) .AND. & - ((amb_info%torsion_b(k)) == (name_atm_b)) .AND. & - ((amb_info%torsion_c(k)) == (name_atm_c)) .AND. & - ((amb_info%torsion_d(k)) == ("X"))) .OR. & - (((amb_info%torsion_a(k)) == ("X")) .AND. & - ((amb_info%torsion_b(k)) == (name_atm_c)) .AND. & - ((amb_info%torsion_c(k)) == (name_atm_b)) .AND. & - ((amb_info%torsion_d(k)) == ("X")))) THEN - 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%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%phi0(imul) = amb_info%torsion_phi0(k) - torsion_list(j)%torsion_kind%nmul = imul - found = .TRUE. - END IF END DO + END IF + END IF ! always have the input param last to overwrite all the other ones @@ -3455,5 +3454,41 @@ CONTAINS END SUBROUTINE store_FF_missing_par +! ************************************************************************************************** +!> \brief Search sorted 2d array of integers for a first occurence of value `val` in row `row` +!> \param array 2d array of integers +!> \param val value to search +!> \param row row to search, default = 1 +!> \return column index if `val` is found in the row `row` of `array`; zero otherwise +! ************************************************************************************************** + FUNCTION bsearch_leftmost_2d(array, val, row) RESULT(res) + INTEGER, INTENT(IN) :: array(:, :), val + INTEGER, INTENT(IN), OPTIONAL :: row + INTEGER :: res + + INTEGER :: left, locRow, mid, right + + locRow = 1 + IF (PRESENT(row)) locRow = row + + left = 1 + right = UBOUND(array, dim=2) + + DO WHILE (left < right) + mid = (left + right)/2 + IF (array(locRow, mid) < val) THEN + left = mid + 1 + ELSE + right = mid + END IF + END DO + + res = left + + ! Not found: + IF (array(locRow, res) /= val) res = 0 + + END FUNCTION bsearch_leftmost_2d + END MODULE force_fields_all diff --git a/src/force_fields_ext.F b/src/force_fields_ext.F index d6fff9160d..4de35d3273 100644 --- a/src/force_fields_ext.F +++ b/src/force_fields_ext.F @@ -762,6 +762,13 @@ CONTAINS amb_info%torsion_k(i) = cp_unit_to_cp2k(amb_info%torsion_k(i), "kcalmol") END DO + DO i = 1, SIZE(amb_info%raw_torsion_k) + + ! Do some units conversion into internal atomic units + amb_info%raw_torsion_phi0(i) = cp_unit_to_cp2k(amb_info%raw_torsion_phi0(i), "rad") + amb_info%raw_torsion_k(i) = cp_unit_to_cp2k(amb_info%raw_torsion_k(i), "kcalmol") + END DO + !----------------------------------------------------------------------------- ! 4. Converts all the Nonbonded info from the param file here !----------------------------------------------------------------------------- diff --git a/src/topology_amber.F b/src/topology_amber.F index c159720f8a..91ebe18237 100644 --- a/src/topology_amber.F +++ b/src/topology_amber.F @@ -556,12 +556,24 @@ CONTAINS CASE ("DIHEDRAL_FORCE_CONSTANT") IF (.NOT. do_forcefield) CYCLE CALL rd_amber_section(parser, section, pk, nptra) + IF (nptra <= 0) CYCLE + ! Save raw values + IF (ASSOCIATED(amb_info%raw_torsion_k)) DEALLOCATE (amb_info%raw_torsion_k) + ALLOCATE (amb_info%raw_torsion_k(nptra), source=pk) CASE ("DIHEDRAL_PERIODICITY") IF (.NOT. do_forcefield) CYCLE CALL rd_amber_section(parser, section, pn, nptra) + IF (nptra <= 0) CYCLE + ! Save raw values + IF (ASSOCIATED(amb_info%raw_torsion_m)) DEALLOCATE (amb_info%raw_torsion_m) + ALLOCATE (amb_info%raw_torsion_m(nptra), source=pn) CASE ("DIHEDRAL_PHASE") IF (.NOT. do_forcefield) CYCLE CALL rd_amber_section(parser, section, phase, nptra) + IF (nptra <= 0) CYCLE + ! Save raw values + IF (ASSOCIATED(amb_info%raw_torsion_phi0)) DEALLOCATE (amb_info%raw_torsion_phi0) + ALLOCATE (amb_info%raw_torsion_phi0(nptra), source=phase) CASE ("LENNARD_JONES_ACOEF") IF (.NOT. do_forcefield) CYCLE CALL rd_amber_section(parser, section, cn1, ntypes*(ntypes + 1)/2) @@ -620,6 +632,25 @@ CONTAINS ! Just Ignore other sections... END SELECT END DO + ! Save raw torsion info: atom indices and dihedral index + IF (do_forcefield .AND. (nphih + nphia > 0)) THEN + IF (ASSOCIATED(amb_info%raw_torsion_id)) DEALLOCATE (amb_info%raw_torsion_id) + ALLOCATE (amb_info%raw_torsion_id(5, nphih + nphia)) + DO i = 1, nphih + amb_info%raw_torsion_id(1, i) = iph(i) + amb_info%raw_torsion_id(2, i) = jph(i) + amb_info%raw_torsion_id(3, i) = kph(i) + amb_info%raw_torsion_id(4, i) = lph(i) + amb_info%raw_torsion_id(5, i) = icph(i) + END DO + DO i = 1, nphia + amb_info%raw_torsion_id(1, nphih + i) = ip(i) + amb_info%raw_torsion_id(2, nphih + i) = jp(i) + amb_info%raw_torsion_id(3, nphih + i) = kp(i) + amb_info%raw_torsion_id(4, nphih + i) = lp(i) + amb_info%raw_torsion_id(5, nphih + i) = icp(i) + END DO + END IF END IF ! Extracts connectivity info from the AMBER topology file @@ -831,6 +862,7 @@ CONTAINS ! Force Fields informations related to torsions ! 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) CALL reallocate(amb_info%torsion_c, 1, buffer_size) @@ -858,6 +890,13 @@ CONTAINS CALL reallocate(amb_info%torsion_m, 1, nsize) CALL reallocate(amb_info%torsion_phi0, 1, nsize) + ! Sort dihedral metadata for faster lookup + IF (nphih + nphia /= 0) THEN + ALLOCATE (iwork(nphih + nphia)) + CALL sort(amb_info%raw_torsion_id, 1, nphih + nphia, 1, 5, iwork) + DEALLOCATE (iwork) + END IF + ! ---------------------------------------------------------- ! Post process of LJ parameters ! ---------------------------------------------------------- @@ -1585,7 +1624,7 @@ CONTAINS label_b(iphi) = work_label(2, 1) label_c(iphi) = work_label(3, 1) label_d(iphi) = work_label(4, 1) - k(iphi) = pk(icp(iwork(1)))*0.5_dp + k(iphi) = pk(icp(iwork(1))) m(iphi) = NINT(pn(icp(iwork(1)))) IF (m(iphi) - pn(icp(iwork(1))) .GT. EPSILON(1.0_dp)) THEN ! non integer torsions not supported @@ -1619,7 +1658,7 @@ CONTAINS label_b(iphi) = work_label(2, i) label_c(iphi) = work_label(3, i) label_d(iphi) = work_label(4, i) - k(iphi) = pk(icp(iwork(i)))*0.5_dp + k(iphi) = pk(icp(iwork(i))) m(iphi) = NINT(pn(icp(iwork(i)))) IF (m(iphi) - pn(icp(iwork(i))) .GT. EPSILON(1.0_dp)) THEN ! non integer torsions not supported diff --git a/tests/Fist/regtest-13/TEST_FILES b/tests/Fist/regtest-13/TEST_FILES index abaad385ce..5b676d372a 100644 --- a/tests/Fist/regtest-13/TEST_FILES +++ b/tests/Fist/regtest-13/TEST_FILES @@ -1,5 +1,5 @@ # multiple_unit_cell with AMBER connectivity -pc-222.inp 2 1.0E-14 -0.105177118522E+02 +pc-222.inp 2 1.0E-14 -0.176219641248E+02 check_ex_14.inp 11 1.0E-14 0.005020140884419 test_ex_14.inp 11 1.0E-14 0.005020140884419 si_muc_cell_opt.inp 7 1.0E-14 -2.7222286209 diff --git a/tests/Fist/regtest-14/TEST_FILES b/tests/Fist/regtest-14/TEST_FILES index c802831fac..2f67d941f8 100644 --- a/tests/Fist/regtest-14/TEST_FILES +++ b/tests/Fist/regtest-14/TEST_FILES @@ -7,7 +7,7 @@ argon-gle_r.inp 2 2.0E-11 - #new FF water_mm3.inp 11 1.0E-14 0.003534941656640 #mol_set + amber -lamol.inp 2 1.0E-14 0.734695738723E+00 +lamol.inp 2 1.0E-14 0.808816509404E+00 #impropers read from input file section imp_test_11.inp 11 1.0E-14 0.804574647188983 imp_test_12.inp 11 1.0E-14 0.804574647188983 diff --git a/tests/Fist/regtest-3/TEST_FILES b/tests/Fist/regtest-3/TEST_FILES index 03e7caa54b..e12e62d836 100644 --- a/tests/Fist/regtest-3/TEST_FILES +++ b/tests/Fist/regtest-3/TEST_FILES @@ -27,7 +27,7 @@ argon.inp 25 1.0E-14 #amber files ace_ala_nme.inp 7 1.0E-14 -0.1185300854 ace_ala_nme-ambconn.inp 7 1.0E-14 -0.1185300850 -ace_ala_nme-amber.inp 7 1.0E-14 -0.1257394724 +ace_ala_nme-amber.inp 7 1.0E-14 -0.1190186541 #EAM alloys agpt.inp 2 1.0E-14 -0.292077327941E+01 #BMHFTD diff --git a/tests/QMMM/SE/regtest-force-mixing/TEST_FILES b/tests/QMMM/SE/regtest-force-mixing/TEST_FILES index 7bfc501fbb..3926fd0843 100644 --- a/tests/QMMM/SE/regtest-force-mixing/TEST_FILES +++ b/tests/QMMM/SE/regtest-force-mixing/TEST_FILES @@ -4,6 +4,6 @@ # 1 compares the last total energy in the file # for details see cp2k/tools/do_regtest # QM/MM -Lysozyme_small_NVT.inp 2 1.0E-14 -0.664207876421E+03 +Lysozyme_small_NVT.inp 2 1.0E-14 -0.663397284084E+03 tyrosine_NVT.inp 2 1.0E-14 -481.71605988300001 #EOF