From 331d841867d89dedc83ab43480a68196364afcfd Mon Sep 17 00:00:00 2001 From: Vladimir Mironov Date: Fri, 11 Sep 2020 15:30:15 +0300 Subject: [PATCH] Use raw Amber parameters for dihedrals Previously, Amber dihedral parameters were assigned in two steps: 1. Find all unique dihedral parameters in the topology. Uniqueness is defined by combination of kinds of involved atoms 2. Walk along the full dihedral list and for every dihedral pick matching set of parameters based on the kinds of corresponding atoms This patch change this behavior to use dihedral parameter assignment found in Amber parameter+topology file. It uses absolute atomic indices instead of kinds and does not lead to disambiguation. --- src/force_field_types.F | 20 +++ src/force_fields_all.F | 115 ++++++++++++------ src/force_fields_ext.F | 7 ++ src/topology_amber.F | 43 ++++++- tests/Fist/regtest-13/TEST_FILES | 2 +- tests/Fist/regtest-14/TEST_FILES | 2 +- tests/Fist/regtest-3/TEST_FILES | 2 +- tests/QMMM/SE/regtest-force-mixing/TEST_FILES | 2 +- 8 files changed, 147 insertions(+), 46 deletions(-) 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