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

@@ -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 @@ SUBROUTINE init_amber_info(amb_info)
       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 @@ SUBROUTINE deallocate_amb_info(amb_info)
       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
 

src/force_fields_all.F

@@ -439,14 +439,19 @@ SUBROUTINE force_field_unique_tors(particle_set, &
                                                             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
       TYPE(torsion_kind_type), DIMENSION(:), POINTER     :: torsion_kind_set
       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 @@ SUBROUTINE force_field_unique_tors(particle_set, &
                           ((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 @@ SUBROUTINE force_field_pack_tors(particle_set, molecule_kind_set, molecule_set,
       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 @@ SUBROUTINE force_field_pack_tors(particle_set, molecule_kind_set, molecule_set,
                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 @@ SUBROUTINE store_FF_missing_par(atm1, atm2, atm3, atm4, type_name, fatal, array)
 
    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
 

src/force_fields_ext.F

@@ -762,6 +762,13 @@ SUBROUTINE read_force_field_amber(ff_type, para_env, mm_section, particle_set)
          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
       !-----------------------------------------------------------------------------

src/topology_amber.F

@@ -556,12 +556,24 @@ SUBROUTINE rdparm_amber_8(filename, output_unit, para_env, do_connectivity, &
             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 @@ SUBROUTINE rdparm_amber_8(filename, output_unit, para_env, do_connectivity, &
                ! 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 @@ SUBROUTINE rdparm_amber_8(filename, output_unit, para_env, do_connectivity, &
          ! 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 @@ SUBROUTINE rdparm_amber_8(filename, output_unit, para_env, do_connectivity, &
          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 @@ SUBROUTINE post_process_torsions_info(label_a, label_b, label_c, label_d, k, &
          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 @@ SUBROUTINE post_process_torsions_info(label_a, label_b, label_c, label_d, k, &
                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

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

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

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

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