mirror of
https://github.com/cp2k/cp2k.git
synced 2026-07-29 06:35:28 -04:00
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.
This commit is contained in:
Vladimir Mironov
GitHub
parent
538497f832
commit
331d841867
8 changed files with 147 additions and 46 deletions
|
|
@ -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
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
!-----------------------------------------------------------------------------
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue