diff --git a/src/force_fields_all.F b/src/force_fields_all.F index d2c0800695..dd3f8e4295 100644 --- a/src/force_fields_all.F +++ b/src/force_fields_all.F @@ -90,6 +90,8 @@ MODULE force_fields_all CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'force_fields_all' PRIVATE + LOGICAL, PARAMETER :: debug_this_module = .FALSE. + PUBLIC :: force_field_unique_bond, & force_field_unique_bend, & force_field_unique_ub, & @@ -121,14 +123,15 @@ CONTAINS !> \param molecule_kind_set ... !> \param molecule_set ... !> \param ff_type ... +!> \param iw ... ! ************************************************************************************************** - SUBROUTINE force_field_unique_bond(particle_set, & - molecule_kind_set, molecule_set, ff_type) + SUBROUTINE force_field_unique_bond(particle_set, molecule_kind_set, molecule_set, ff_type, iw) TYPE(particle_type), DIMENSION(:), POINTER :: particle_set TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set TYPE(force_field_type), INTENT(INOUT) :: ff_type + INTEGER, INTENT(IN) :: iw CHARACTER(len=*), PARAMETER :: routineN = 'force_field_unique_bond' @@ -146,6 +149,12 @@ CONTAINS TYPE(molecule_type), POINTER :: molecule CALL timeset(routineN, handle2) + + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(/,T2,A)") & + "FORCEFIELD| Checking for unique bond terms" + END IF + DO i = 1, SIZE(molecule_kind_set) molecule_kind => molecule_kind_set(i) CALL get_molecule_kind(molecule_kind=molecule_kind, & @@ -217,14 +226,15 @@ CONTAINS !> \param molecule_kind_set ... !> \param molecule_set ... !> \param ff_type ... +!> \param iw ... ! ************************************************************************************************** - SUBROUTINE force_field_unique_bend(particle_set, & - molecule_kind_set, molecule_set, ff_type) + SUBROUTINE force_field_unique_bend(particle_set, molecule_kind_set, molecule_set, ff_type, iw) TYPE(particle_type), DIMENSION(:), POINTER :: particle_set TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set TYPE(force_field_type), INTENT(INOUT) :: ff_type + INTEGER, INTENT(IN) :: iw CHARACTER(len=*), PARAMETER :: routineN = 'force_field_unique_bend' @@ -242,6 +252,12 @@ CONTAINS TYPE(molecule_type), POINTER :: molecule CALL timeset(routineN, handle2) + + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(/,T2,A)") & + "FORCEFIELD| Checking for unique bend terms" + END IF + DO i = 1, SIZE(molecule_kind_set) molecule_kind => molecule_kind_set(i) CALL get_molecule_kind(molecule_kind=molecule_kind, & @@ -313,6 +329,7 @@ CONTAINS DEALLOCATE (map_bend_kind) END IF END DO + CALL timestop(handle2) END SUBROUTINE force_field_unique_bend @@ -322,13 +339,14 @@ CONTAINS !> \param particle_set ... !> \param molecule_kind_set ... !> \param molecule_set ... +!> \param iw ... ! ************************************************************************************************** - SUBROUTINE force_field_unique_ub(particle_set, & - molecule_kind_set, molecule_set) + SUBROUTINE force_field_unique_ub(particle_set, molecule_kind_set, molecule_set, iw) TYPE(particle_type), DIMENSION(:), POINTER :: particle_set TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set + INTEGER, INTENT(IN) :: iw CHARACTER(len=*), PARAMETER :: routineN = 'force_field_unique_ub' @@ -346,6 +364,12 @@ CONTAINS TYPE(ub_type), DIMENSION(:), POINTER :: ub_list CALL timeset(routineN, handle2) + + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(/,T2,A)") & + "FORCEFIELD| Checking for unique Urey-Bradley terms" + END IF + DO i = 1, SIZE(molecule_kind_set) molecule_kind => molecule_kind_set(i) CALL get_molecule_kind(molecule_kind=molecule_kind, & @@ -419,14 +443,15 @@ CONTAINS !> \param molecule_kind_set ... !> \param molecule_set ... !> \param ff_type ... +!> \param iw ... ! ************************************************************************************************** - SUBROUTINE force_field_unique_tors(particle_set, & - molecule_kind_set, molecule_set, ff_type) + SUBROUTINE force_field_unique_tors(particle_set, molecule_kind_set, molecule_set, ff_type, iw) TYPE(particle_type), DIMENSION(:), POINTER :: particle_set TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set TYPE(force_field_type), INTENT(INOUT) :: ff_type + INTEGER, INTENT(IN) :: iw CHARACTER(len=*), PARAMETER :: routineN = 'force_field_unique_tors' @@ -447,6 +472,11 @@ CONTAINS CALL timeset(routineN, handle2) + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(/,T2,A)") & + "FORCEFIELD| Checking for unique torsion terms" + END IF + ! 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) @@ -533,6 +563,7 @@ CONTAINS DEALLOCATE (map_torsion_kind) END IF END DO + CALL timestop(handle2) END SUBROUTINE force_field_unique_tors @@ -543,14 +574,15 @@ CONTAINS !> \param molecule_kind_set ... !> \param molecule_set ... !> \param ff_type ... +!> \param iw ... ! ************************************************************************************************** - SUBROUTINE force_field_unique_impr(particle_set, & - molecule_kind_set, molecule_set, ff_type) + SUBROUTINE force_field_unique_impr(particle_set, molecule_kind_set, molecule_set, ff_type, iw) TYPE(particle_type), DIMENSION(:), POINTER :: particle_set TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set TYPE(force_field_type), INTENT(INOUT) :: ff_type + INTEGER, INTENT(IN) :: iw CHARACTER(len=*), PARAMETER :: routineN = 'force_field_unique_impr' @@ -570,6 +602,12 @@ CONTAINS TYPE(molecule_type), POINTER :: molecule CALL timeset(routineN, handle2) + + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(/,T2,A)") & + "FORCEFIELD| Checking for unique improper terms" + END IF + DO i = 1, SIZE(molecule_kind_set) molecule_kind => molecule_kind_set(i) CALL get_molecule_kind(molecule_kind=molecule_kind, & @@ -577,7 +615,9 @@ CONTAINS natom=natom, & nimpr=nimpr, impr_list=impr_list) molecule => molecule_set(molecule_list(1)) + CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last) + IF (nimpr > 0) THEN ALLOCATE (map_impr_kind(nimpr)) counter = 0 @@ -663,14 +703,15 @@ CONTAINS !> \param molecule_kind_set ... !> \param molecule_set ... !> \param ff_type ... +!> \param iw ... ! ************************************************************************************************** - SUBROUTINE force_field_unique_opbend(particle_set, & - molecule_kind_set, molecule_set, ff_type) + SUBROUTINE force_field_unique_opbend(particle_set, molecule_kind_set, molecule_set, ff_type, iw) TYPE(particle_type), DIMENSION(:), POINTER :: particle_set TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set TYPE(force_field_type), INTENT(INOUT) :: ff_type + INTEGER, INTENT(IN) :: iw CHARACTER(len=*), PARAMETER :: routineN = 'force_field_unique_opbend' @@ -690,6 +731,12 @@ CONTAINS TYPE(opbend_type), DIMENSION(:), POINTER :: opbend_list CALL timeset(routineN, handle2) + + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(/,T2,A)") & + "FORCEFIELD| Checking for unique out-of-plane bend terms" + END IF + DO i = 1, SIZE(molecule_kind_set) molecule_kind => molecule_kind_set(i) CALL get_molecule_kind(molecule_kind=molecule_kind, & @@ -786,9 +833,10 @@ CONTAINS !> \param inp_info ... !> \param gro_info ... !> \param amb_info ... +!> \param iw ... ! ************************************************************************************************** - SUBROUTINE force_field_pack_bond(particle_set, molecule_kind_set, molecule_set, & - fatal, Ainfo, chm_info, inp_info, gro_info, amb_info) + SUBROUTINE force_field_pack_bond(particle_set, molecule_kind_set, molecule_set, fatal, Ainfo, & + chm_info, inp_info, gro_info, amb_info, iw) TYPE(particle_type), DIMENSION(:), POINTER :: particle_set TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set @@ -800,6 +848,7 @@ CONTAINS TYPE(input_info_type), POINTER :: inp_info TYPE(gromos_info_type), POINTER :: gro_info TYPE(amber_info_type), POINTER :: amb_info + INTEGER, INTENT(IN) :: iw CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_bond' @@ -815,6 +864,11 @@ CONTAINS CALL timeset(routineN, handle2) + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(/,T2,A)") & + "FORCEFIELD| Checking for bond terms" + END IF + DO i = 1, SIZE(molecule_kind_set) molecule_kind => molecule_kind_set(i) CALL get_molecule_kind(molecule_kind=molecule_kind, & @@ -917,9 +971,12 @@ CONTAINS bond_list(j)%bond_kind%id_type = do_ff_undef END IF END DO + CALL set_molecule_kind(molecule_kind=molecule_kind, & bond_list=bond_list) + END DO + CALL timestop(handle2) END SUBROUTINE force_field_pack_bond @@ -935,9 +992,10 @@ CONTAINS !> \param inp_info ... !> \param gro_info ... !> \param amb_info ... +!> \param iw ... ! ************************************************************************************************** - SUBROUTINE force_field_pack_bend(particle_set, molecule_kind_set, molecule_set, & - fatal, Ainfo, chm_info, inp_info, gro_info, amb_info) + SUBROUTINE force_field_pack_bend(particle_set, molecule_kind_set, molecule_set, fatal, Ainfo, & + chm_info, inp_info, gro_info, amb_info, iw) TYPE(particle_type), DIMENSION(:), POINTER :: particle_set TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set @@ -949,6 +1007,7 @@ CONTAINS TYPE(input_info_type), POINTER :: inp_info TYPE(gromos_info_type), POINTER :: gro_info TYPE(amber_info_type), POINTER :: amb_info + INTEGER, INTENT(IN) :: iw CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_bend' @@ -964,6 +1023,11 @@ CONTAINS CALL timeset(routineN, handle2) + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(/,T2,A)") & + "FORCEFIELD| Checking for bend terms" + END IF + DO i = 1, SIZE(molecule_kind_set) molecule_kind => molecule_kind_set(i) CALL get_molecule_kind(molecule_kind=molecule_kind, & @@ -1137,6 +1201,12 @@ CONTAINS TYPE(ub_type), DIMENSION(:), POINTER :: ub_list CALL timeset(routineN, handle2) + + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(/,T2,A)") & + "FORCEFIELD| Checking for Urey-Bradley (UB) terms" + END IF + DO i = 1, SIZE(molecule_kind_set) molecule_kind => molecule_kind_set(i) CALL get_molecule_kind(molecule_kind=molecule_kind, & @@ -1164,10 +1234,10 @@ CONTAINS CALL uppercase(name_atm_b) CALL uppercase(name_atm_c) - ! loop over params from GROMOS + ! Loop over params from GROMOS ! ikuo - None that I know... - ! loop over params from CHARMM + ! Loop over params from CHARMM IF (ASSOCIATED(chm_info%ub_a)) THEN DO k = 1, SIZE(chm_info%ub_a) IF ((((chm_info%ub_a(k)) == (name_atm_a)) .AND. & @@ -1179,8 +1249,11 @@ CONTAINS ub_list(j)%ub_kind%id_type = do_ff_charmm ub_list(j)%ub_kind%k(1) = chm_info%ub_k(k) ub_list(j)%ub_kind%r0 = chm_info%ub_r0(k) - IF (iw > 0) WRITE (iw, *) " Found UB ", TRIM(name_atm_a), " ", & - TRIM(name_atm_b), " ", TRIM(name_atm_c) + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(T2,A)") & + "FORCEFIELD| Found Urey-Bradley term (CHARMM) for the atomic kinds "// & + TRIM(name_atm_a)//", "//TRIM(name_atm_b)//" and "//TRIM(name_atm_c) + END IF CALL issue_duplications(found, "Urey-Bradley", name_atm_a, & name_atm_b, name_atm_c) found = .TRUE. @@ -1189,10 +1262,10 @@ CONTAINS END DO END IF - ! loop over params from AMBER + ! Loop over params from AMBER ! teo - None that I know... - ! always have the input param last to overwrite all the other ones + ! Always have the input param last to overwrite all the other ones IF (ASSOCIATED(inp_info%ub_a)) THEN DO k = 1, SIZE(inp_info%ub_a) IF ((((inp_info%ub_a(k)) == (name_atm_a)) .AND. & @@ -1204,6 +1277,11 @@ CONTAINS ub_list(j)%ub_kind%id_type = inp_info%ub_kind(k) ub_list(j)%ub_kind%k(:) = inp_info%ub_k(:, k) ub_list(j)%ub_kind%r0 = inp_info%ub_r0(k) + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(T2,A)") & + "FORCEFIELD| Found Urey-Bradley term (input) for the atomic kinds "// & + TRIM(name_atm_a)//", "//TRIM(name_atm_b)//" and "//TRIM(name_atm_c) + END IF CALL issue_duplications(found, "Urey-Bradley", name_atm_a, & name_atm_b, name_atm_c) found = .TRUE. @@ -1230,9 +1308,12 @@ CONTAINS ub_list(j)%ub_kind%id_type = do_ff_undef END IF END DO + CALL set_molecule_kind(molecule_kind=molecule_kind, & ub_list=ub_list) + END DO + CALL timestop(handle2) END SUBROUTINE force_field_pack_ub @@ -1265,8 +1346,8 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_tors' - CHARACTER(LEN=default_string_length) :: ldum, name_atm_a, name_atm_b, & - name_atm_c, name_atm_d + CHARACTER(LEN=default_string_length) :: ldum, molecule_name, 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, k_end, & k_start, last, natom, ntorsion, & @@ -1281,14 +1362,23 @@ CONTAINS CALL timeset(routineN, handle2) + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(/,T2,A)") & + "FORCEFIELD| Checking for torsion terms" + END IF + DO i = 1, SIZE(molecule_kind_set) molecule_kind => molecule_kind_set(i) CALL get_molecule_kind(molecule_kind=molecule_kind, & molecule_list=molecule_list, & + name=molecule_name, & natom=natom, & - ntorsion=ntorsion, torsion_list=torsion_list) + ntorsion=ntorsion, & + torsion_list=torsion_list) molecule => molecule_set(molecule_list(1)) - CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last) + CALL get_molecule(molecule=molecule, & + first_atom=first, & + last_atom=last) DO j = 1, ntorsion IF (torsion_list(j)%torsion_kind%id_type == do_ff_undef) THEN atm_a = torsion_list(j)%a @@ -1314,7 +1404,7 @@ CONTAINS CALL uppercase(name_atm_c) CALL uppercase(name_atm_d) - ! loop over params from GROMOS + ! Loop over params from GROMOS IF (ASSOCIATED(gro_info%torsion_k)) THEN k = SIZE(gro_info%torsion_k) itype = torsion_list(j)%itype @@ -1341,7 +1431,7 @@ CONTAINS imul = torsion_list(j)%torsion_kind%nmul END IF - ! loop over params from CHARMM + ! Loop over params from CHARMM IF (ASSOCIATED(chm_info%torsion_a)) THEN DO k = 1, SIZE(chm_info%torsion_a) IF ((((chm_info%torsion_a(k)) == (name_atm_a)) .AND. & @@ -1390,7 +1480,7 @@ CONTAINS END IF END IF - ! loop over params from AMBER + ! 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 @@ -1430,7 +1520,7 @@ CONTAINS END IF - ! always have the input param last to overwrite all the other ones + ! Always have the input param last to overwrite all the other ones IF (ASSOCIATED(inp_info%torsion_a)) THEN DO k = 1, SIZE(inp_info%torsion_a) IF ((((inp_info%torsion_a(k)) == (name_atm_a)) .AND. & @@ -1455,7 +1545,28 @@ CONTAINS END DO END IF - IF (.NOT. found) THEN + IF (found) THEN + ldum = cp_to_string(imul) + IF (iw > 0) THEN + IF (imul < 1) THEN + WRITE (UNIT=iw, FMT="(T2,A)") & + "FORCEFIELD| No torsion term found" + ELSE IF (imul == 1) THEN + WRITE (UNIT=iw, FMT="(T2,A)") & + "FORCEFIELD| Found torsion term for the atomic kinds "// & + TRIM(name_atm_a)//", "//TRIM(name_atm_b)// & + ", "//TRIM(name_atm_c)// & + " and "//TRIM(name_atm_d) + ELSE + WRITE (UNIT=iw, FMT="(T2,A)") & + "FORCEFIELD| Found multiple ("//TRIM(ldum)// & + ") torsion terms for the atomic kinds "// & + TRIM(name_atm_a)//", "//TRIM(name_atm_b)// & + ", "//TRIM(name_atm_c)// & + " and "//TRIM(name_atm_d) + END IF + END IF + ELSE CALL store_FF_missing_par(atm1=TRIM(name_atm_a), & atm2=TRIM(name_atm_b), & atm3=TRIM(name_atm_c), & @@ -1464,35 +1575,32 @@ CONTAINS array=Ainfo) torsion_list(j)%torsion_kind%id_type = do_ff_undef torsion_list(j)%id_type = do_ff_undef - ELSE - ldum = cp_to_string(imul) - IF ((imul /= 1) .AND. (iw > 0)) & - WRITE (iw, '(/,2("UTIL_INFO| ",A,/))') & - "Multiple Torsion declarations: "//TRIM(name_atm_a)// & - ","//TRIM(name_atm_b)//","//TRIM(name_atm_c)//" and "//TRIM(name_atm_d), & - "Number of defined torsions: "//TRIM(ldum)//" ." END IF - ! + ! QM/MM modifications - ! IF (only_qm) THEN - IF (iw > 0) WRITE (iw, *) " Torsion PARAM between QM atoms ", j, " : ", & - TRIM(name_atm_a), " ", & - TRIM(name_atm_b), " ", & - TRIM(name_atm_c), " ", & - TRIM(name_atm_d), " ", & - torsion_list(j)%a, & - torsion_list(j)%b, & - torsion_list(j)%c, & - torsion_list(j)%d + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(T2,A,I0,4(A,I0))") & + "FORCEFIELD| Torsion ", j, " for molecule kind "//TRIM(molecule_name)// & + TRIM(name_atm_a)// & + "-"//TRIM(name_atm_b)//"-"//TRIM(name_atm_c)//"-"// & + TRIM(name_atm_d)//" (", torsion_list(j)%a, ", ", & + torsion_list(j)%b, ", ", torsion_list(j)%c, ", ", & + torsion_list(j)%d + END IF torsion_list(j)%torsion_kind%id_type = do_ff_undef torsion_list(j)%id_type = do_ff_undef END IF + END IF - END DO + + END DO ! torsion + CALL set_molecule_kind(molecule_kind=molecule_kind, & torsion_list=torsion_list) - END DO + + END DO ! molecule kind + CALL timestop(handle2) END SUBROUTINE force_field_pack_tors @@ -1506,9 +1614,10 @@ CONTAINS !> \param chm_info ... !> \param inp_info ... !> \param gro_info ... +!> \param iw ... ! ************************************************************************************************** SUBROUTINE force_field_pack_impr(particle_set, molecule_kind_set, molecule_set, & - Ainfo, chm_info, inp_info, gro_info) + Ainfo, chm_info, inp_info, gro_info, iw) TYPE(particle_type), DIMENSION(:), POINTER :: particle_set TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set @@ -1518,6 +1627,7 @@ CONTAINS TYPE(charmm_info_type), POINTER :: chm_info TYPE(input_info_type), POINTER :: inp_info TYPE(gromos_info_type), POINTER :: gro_info + INTEGER, INTENT(IN) :: iw CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_impr' @@ -1535,14 +1645,23 @@ CONTAINS CALL timeset(routineN, handle2) + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(/,T2,A)") & + "FORCEFIELD| Checking for improper terms" + END IF + DO i = 1, SIZE(molecule_kind_set) + molecule_kind => molecule_kind_set(i) CALL get_molecule_kind(molecule_kind=molecule_kind, & molecule_list=molecule_list, & natom=natom, & - nimpr=nimpr, impr_list=impr_list) + nimpr=nimpr, & + impr_list=impr_list) + molecule => molecule_set(molecule_list(1)) CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last) + DO j = 1, nimpr atm_a = impr_list(j)%a atomic_kind => particle_set(atm_a + first - 1)%atomic_kind @@ -1567,7 +1686,7 @@ CONTAINS CALL uppercase(name_atm_c) CALL uppercase(name_atm_d) - ! loop over params from GROMOS + ! Loop over params from GROMOS IF (ASSOCIATED(gro_info%impr_k)) THEN k = SIZE(gro_info%impr_k) itype = impr_list(j)%itype @@ -1583,7 +1702,7 @@ CONTAINS impr_list(j)%id_type = gro_info%ff_gromos_type END IF - ! loop over params from CHARMM + ! Loop over params from CHARMM IF (ASSOCIATED(chm_info%impr_a)) THEN DO k = 1, SIZE(chm_info%impr_a) IF ((((chm_info%impr_a(k)) == (name_atm_a)) .AND. & @@ -1625,7 +1744,7 @@ CONTAINS END IF END IF - ! loop over params from AMBER not needed since impropers in AMBER + ! Loop over params from AMBER not needed since impropers in AMBER ! are treated like standard torsions ! always have the input param last to overwrite all the other ones @@ -1674,17 +1793,27 @@ CONTAINS impr_list(j)%impr_kind%id_type = do_ff_undef impr_list(j)%id_type = do_ff_undef END IF - ! + ! QM/MM modifications - ! IF (only_qm) THEN + IF (found) THEN + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(T2,A)") & + "FORCEFIELD| Found improper term for "//TRIM(name_atm_a)// & + "-"//TRIM(name_atm_b)//"-"//TRIM(name_atm_c)//"-"// & + TRIM(name_atm_d) + END IF + END IF impr_list(j)%impr_kind%id_type = do_ff_undef impr_list(j)%id_type = do_ff_undef END IF END DO + CALL set_molecule_kind(molecule_kind=molecule_kind, impr_list=impr_list) + END DO + CALL timestop(handle2) END SUBROUTINE force_field_pack_impr @@ -1698,10 +1827,11 @@ CONTAINS !> \param molecule_set ... !> \param Ainfo ... !> \param inp_info ... +!> \param iw ... !> \author Louis Vanduyfhuys ! ************************************************************************************************** - SUBROUTINE force_field_pack_opbend(particle_set, molecule_kind_set, molecule_set, & - Ainfo, inp_info) + SUBROUTINE force_field_pack_opbend(particle_set, molecule_kind_set, molecule_set, Ainfo, & + inp_info, iw) TYPE(particle_type), DIMENSION(:), POINTER :: particle_set TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set @@ -1709,6 +1839,7 @@ CONTAINS CHARACTER(LEN=default_string_length), & DIMENSION(:), POINTER :: Ainfo TYPE(input_info_type), POINTER :: inp_info + INTEGER, INTENT(IN) :: iw CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_opbend' @@ -1725,6 +1856,11 @@ CONTAINS CALL timeset(routineN, handle2) + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(/,T2,A)") & + "FORCEFIELD| Checking for out-of-plane bend terms" + END IF + DO i = 1, SIZE(molecule_kind_set) molecule_kind => molecule_kind_set(i) CALL get_molecule_kind(molecule_kind=molecule_kind, & @@ -1812,8 +1948,11 @@ CONTAINS END IF END DO + CALL set_molecule_kind(molecule_kind=molecule_kind, opbend_list=opbend_list) + END DO + CALL timestop(handle2) END SUBROUTINE force_field_pack_opbend @@ -1839,7 +1978,7 @@ CONTAINS LOGICAL :: my_qmmm TYPE(qmmm_env_mm_type), POINTER :: qmmm_env TYPE(input_info_type), POINTER :: inp_info - INTEGER :: iw4 + INTEGER, INTENT(IN) :: iw4 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_charges' @@ -1854,12 +1993,13 @@ CONTAINS TYPE(val_type), POINTER :: val CALL timeset(routineN, handle) + charge_tot = 0.0_dp NULLIFY (list) ! Not implemented for core-shell IF (ASSOCIATED(inp_info%shell_list)) THEN - CPABORT("Array of charges not implemented for CORE-SHELL model!") + CPABORT("Array of charges is not implemented for the core-shell model") END IF ! Allocate array to particle_set size @@ -1918,16 +2058,13 @@ CONTAINS END IF END IF END IF - ! + ! QM/MM modifications - ! IF (only_qm .AND. my_qmmm) THEN IF (qmmm_env%qmmm_coupl_type /= do_qmmm_none) THEN scale_factor = 0.0_dp IF (is_link_atom) THEN - ! ! Find the scaling factor... - ! DO ilink = 1, SIZE(qmmm_env%mm_link_atoms) IF (iatom == qmmm_env%mm_link_atoms(ilink)) EXIT END DO @@ -1938,17 +2075,22 @@ CONTAINS END IF END IF END DO + ! Sum up total charges for IO charge_tot = SUM(charges) - ! Print Total Charge of the system + + ! Print total charge of the system IF (iw4 > 0) THEN - WRITE (iw4, FMT='(T2,"CHARGE_INFO| Total Charge of the Classical System: ",T69,F12.6)') charge_tot + WRITE (UNIT=iw4, FMT="(/,T2,A,T61,F20.10)") & + "FORCEFIELD| Total charge of the classical system: ", charge_tot END IF + CALL timestop(handle) + END SUBROUTINE force_field_pack_charges ! ************************************************************************************************** -!> \brief Set up atomic_kind_set()%fist_potentail%[qeff] +!> \brief Set up atomic_kind_set()%fist_potential%[qeff] !> and shell potential parameters !> \param atomic_kind_set ... !> \param qmmm_env ... @@ -1964,11 +2106,11 @@ CONTAINS TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(qmmm_env_mm_type), POINTER :: qmmm_env - LOGICAL :: fatal - INTEGER :: iw, iw4 + LOGICAL, INTENT(INOUT) :: fatal + INTEGER, INTENT(IN) :: iw, iw4 CHARACTER(LEN=default_string_length), & DIMENSION(:), POINTER :: Ainfo - LOGICAL :: my_qmmm + LOGICAL, INTENT(IN) :: my_qmmm TYPE(input_info_type), POINTER :: inp_info CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_charge' @@ -1984,7 +2126,9 @@ CONTAINS TYPE(fist_potential_type), POINTER :: fist_potential CALL timeset(routineN, handle) + charge_tot = 0.0_dp + DO i = 1, SIZE(atomic_kind_set) atomic_kind => atomic_kind_set(i) CALL get_atomic_kind(atomic_kind=atomic_kind, & @@ -2001,8 +2145,15 @@ CONTAINS ! Always have the input param last to overwrite all the other ones IF (ASSOCIATED(inp_info%charge_atm)) THEN + IF (iw > 0) WRITE (UNIT=iw, FMT="(A)") "" DO j = 1, SIZE(inp_info%charge_atm) - IF (iw > 0) WRITE (iw, *) "Charge Checking ::", TRIM(inp_info%charge_atm(j)), atmname + IF (debug_this_module) THEN + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(T2,A)") & + "Checking charges for the atomic kinds "// & + TRIM(inp_info%charge_atm(j))//" and "//TRIM(atmname) + END IF + END IF IF ((inp_info%charge_atm(j)) == atmname) THEN charge = inp_info%charge(j) CALL issue_duplications(found, "Charge", atmname) @@ -2097,11 +2248,15 @@ CONTAINS END IF END IF END DO - ! Print Total Charge of the system + + ! Print total charge of the system IF (iw4 > 0) THEN - WRITE (iw4, FMT='(T2,"CHARGE_INFO| Total Charge of the Classical System: ",T69,F12.6)') charge_tot + WRITE (UNIT=iw4, FMT="(/,T2,A,T61,F20.10)") & + "FORCEFIELD| Total charge of the classical system: ", charge_tot END IF + CALL timestop(handle) + END SUBROUTINE force_field_pack_charge ! ************************************************************************************************** @@ -2142,13 +2297,17 @@ CONTAINS ! Try to find a matching KIND section in the SUBSYS section and read the ! MM_RADIUS field if it is present. In case the kind section is never ! encountered, the mm_radius remains zero. + IF (iw > 0) WRITE (UNIT=iw, FMT="(A)") "" mm_radius = 0.0_dp DO i_rep = 1, n_rep CALL section_vals_val_get(kind_section, "_SECTION_PARAMETERS_", & c_val=inp_kind_name, i_rep_section=i_rep) CALL uppercase(inp_kind_name) - IF (iw > 0) WRITE (iw, *) "Matching kinds for MM_RADIUS :: '", & - TRIM(kind_name), "' with '", TRIM(inp_kind_name), "'" + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(T2,A)") & + "FORCEFIELD| Matching atomic kinds "//TRIM(kind_name)// & + " and "//TRIM(inp_kind_name)//" for MM_RADIUS" + END IF IF (TRIM(kind_name) == TRIM(inp_kind_name)) THEN CALL section_vals_val_get(kind_section, i_rep_section=i_rep, & keyword_name="MM_RADIUS", r_val=mm_radius) @@ -2156,10 +2315,11 @@ CONTAINS found = .TRUE. END IF END DO - CALL set_potential(potential=fist_potential, mm_radius=mm_radius) END DO + CALL timestop(handle) + END SUBROUTINE force_field_pack_radius ! ************************************************************************************************** @@ -2186,6 +2346,11 @@ CONTAINS CALL timeset(routineN, handle) + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(/,T2,A)") & + "FORCEFIELD| Checking for polarisable forcefield terms" + END IF + DO i = 1, SIZE(atomic_kind_set) atomic_kind => atomic_kind_set(i) CALL get_atomic_kind(atomic_kind=atomic_kind, & @@ -2195,11 +2360,15 @@ CONTAINS CALL uppercase(kind_name) found = .FALSE. + IF (iw > 0) WRITE (UNIT=iw, FMT="(A)") "" ! Always have the input param last to overwrite all the other ones IF (ASSOCIATED(inp_info%apol_atm)) THEN DO j = 1, SIZE(inp_info%apol_atm) - IF (iw > 0) WRITE (iw, *) "Matching kinds for APOL :: '", & - TRIM(kind_name), "' with '", TRIM(inp_info%apol_atm(j)), "'" + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(T2,A)") & + "FORCEFIELD| Matching atomic kinds "//TRIM(kind_name)// & + " and "//TRIM(inp_info%apol_atm(j))//" for APOL" + END IF IF ((inp_info%apol_atm(j)) == kind_name) THEN apol = inp_info%apol(j) CALL issue_duplications(found, "APOL", kind_name) @@ -2210,8 +2379,11 @@ CONTAINS IF (ASSOCIATED(inp_info%cpol_atm)) THEN DO j = 1, SIZE(inp_info%cpol_atm) - IF (iw > 0) WRITE (iw, *) "Matching kinds for CPOL :: '", & - TRIM(kind_name), "' with '", TRIM(inp_info%cpol_atm(j)), "'" + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(T2,A)") & + "FORCEFIELD| Matching atomic kinds "//TRIM(kind_name)// & + " and "//TRIM(inp_info%cpol_atm(j))//" for CPOL" + END IF IF ((inp_info%cpol_atm(j)) == kind_name) THEN cpol = inp_info%cpol(j) CALL issue_duplications(found, "CPOL", kind_name) @@ -2221,8 +2393,11 @@ CONTAINS END IF CALL set_potential(potential=fist_potential, apol=apol, cpol=cpol) + END DO + CALL timestop(handle) + END SUBROUTINE force_field_pack_pol ! ************************************************************************************************** @@ -2231,8 +2406,7 @@ CONTAINS !> \param iw ... !> \param inp_info ... ! ************************************************************************************************** - SUBROUTINE force_field_pack_damp(atomic_kind_set, & - iw, inp_info) + SUBROUTINE force_field_pack_damp(atomic_kind_set, iw, inp_info) TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set INTEGER :: iw @@ -2248,38 +2422,48 @@ CONTAINS TYPE(damping_p_type), POINTER :: damping CALL timeset(routineN, handle2) + + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(/,T2,A)") & + "FORCEFIELD| Checking for damping terms" + END IF + NULLIFY (damping) nkinds = SIZE(atomic_kind_set) DO j = 1, SIZE(atomic_kind_set) + atomic_kind => atomic_kind_set(j) + CALL get_atomic_kind(atomic_kind=atomic_kind, & name=atm_name1) - CALL UPPERCASE(atm_name1) + CALL uppercase(atm_name1) + IF (ASSOCIATED(inp_info%damping_list)) THEN DO i = 1, SIZE(inp_info%damping_list) my_atm_name1 = inp_info%damping_list(i)%atm_name1 my_atm_name2 = inp_info%damping_list(i)%atm_name2 - - IF (iw > 0) WRITE (iw, *) "Damping Checking ::", TRIM(my_atm_name1), & - TRIM(atm_name1) + IF (debug_this_module) THEN + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(T2,A)") & + "FORCEFIELD| Check damping for the atomic kinds "// & + TRIM(my_atm_name1)//" and "//TRIM(atm_name1) + END IF + END IF IF (my_atm_name1 == atm_name1) THEN IF (.NOT. ASSOCIATED(damping)) THEN CALL damping_p_create(damping, nkinds) END IF - found = .FALSE. DO k = 1, SIZE(atomic_kind_set) atomic_kind2 => atomic_kind_set(k) CALL get_atomic_kind(atomic_kind=atomic_kind2, & name=atm_name2) - CALL UPPERCASE(atm_name2) - + CALL uppercase(atm_name2) IF (my_atm_name2 == atm_name2) THEN IF (damping%damp(k)%bij /= HUGE(0.0_dp)) found = .TRUE. CALL issue_duplications(found, "Damping", atm_name1) found = .TRUE. - SELECT CASE (TRIM(inp_info%damping_list(i)%dtype)) CASE ('TANG-TOENNIES') damping%damp(k)%itype = tang_toennies @@ -2290,20 +2474,21 @@ CONTAINS damping%damp(k)%bij = inp_info%damping_list(i)%bij damping%damp(k)%cij = inp_info%damping_list(i)%cij END IF - END DO - IF (.NOT. found) & + IF (.NOT. found) THEN CALL cp_warn(__LOCATION__, & "Atom "//TRIM(my_atm_name2)// & " in damping parameters for atom "//TRIM(my_atm_name1)// & - " not found! ") + " not found.") + END IF END IF END DO - END IF CALL set_atomic_kind(atomic_kind=atomic_kind, damping=damping) + NULLIFY (damping) + END DO CALL timestop(handle2) @@ -2356,6 +2541,7 @@ CONTAINS TYPE(shell_type), DIMENSION(:), POINTER :: shell_list CALL timeset(routineN, handle2) + nshell_tot = 0 n = 0 first_shell = 1 @@ -2367,6 +2553,11 @@ CONTAINS global_section => section_vals_get_subs_vals(root_section, "GLOBAL") CALL section_vals_val_get(global_section, "SAVE_MEM", l_val=save_mem) + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(/,T2,A)") & + "FORCEFIELD| Checking for core-shell terms" + END IF + DO i = 1, SIZE(atomic_kind_set) atomic_kind => atomic_kind_set(i) CALL get_atomic_kind(atomic_kind=atomic_kind, & @@ -2379,7 +2570,13 @@ CONTAINS ! The shell potential can be defined only from input IF (ASSOCIATED(inp_info%shell_list)) THEN DO j = 1, SIZE(inp_info%shell_list) - IF (iw > 0) WRITE (iw, *) "Shell Checking ::", TRIM(inp_info%shell_list(j)%atm_name), atmname + IF (debug_this_module) THEN + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(T2,A)") & + "Checking shells for the atomic kinds "// & + TRIM(inp_info%shell_list(j)%atm_name)//" and "//TRIM(atmname) + END IF + END IF IF ((inp_info%shell_list(j)%atm_name) == atmname) THEN CALL get_atomic_kind(atomic_kind=atomic_kind, & shell=shell, mass=atmmass, natom=natom) @@ -2400,11 +2597,14 @@ CONTAINS CALL set_atomic_kind(atomic_kind=atomic_kind, & shell=shell, shell_active=.TRUE.) END IF - END DO ! j shell kind + END DO ! shell kind END IF ! associated shell_list - END DO ! i atomic kind + END DO ! atomic kind - IF (iw > 0) WRITE (iw, *) "Total number of particles with a shell :: ", nshell_tot + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(/,T2,A,T61,I20)") & + "FORCEFIELD| Total number of particles with a shell:", nshell_tot + END IF ! If shell-model is present: Create particle_set of shells (coord. vel. force) NULLIFY (shell_particle_set) NULLIFY (core_particle_set) @@ -2510,8 +2710,10 @@ CONTAINS n = n + nmol*counter END DO ! i molecule kind END IF + CPASSERT(first_shell - 1 == nshell_tot) CPASSERT(n == nshell_tot) + CALL timestop(handle2) END SUBROUTINE force_field_pack_shell @@ -2558,11 +2760,19 @@ CONTAINS TYPE(atomic_kind_type), POINTER :: atomic_kind TYPE(pair_potential_single_type), POINTER :: pot + CALL timeset(routineN, handle2) + use_qmmm_ff = qmmm_env%use_qmmm_ff NULLIFY (pot) + + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(/,T2,A)") & + "FORCEFIELD| Checking for nonbonded14 terms" + END IF + CALL ewald_env_get(ewald_env, rcut=ewald_rcut) - CALL timeset(routineN, handle2) CALL pair_potential_pp_create(potparm_nonbond14, SIZE(atomic_kind_set)) + DO i = 1, SIZE(atomic_kind_set) atomic_kind => atomic_kind_set(i) CALL get_atomic_kind(atomic_kind=atomic_kind, name=name_atm_a_local) @@ -2611,7 +2821,7 @@ CONTAINS END IF END IF - ! loop over params from CHARMM + ! Loop over params from CHARMM ii = 0 jj = 0 IF (ASSOCIATED(chm_info%nonbond_a_14)) THEN @@ -2668,7 +2878,7 @@ CONTAINS found = .TRUE. END IF - ! loop over params from AMBER + ! Loop over params from AMBER IF (ASSOCIATED(amb_info%nonbond_a)) THEN ii = 0 jj = 0 @@ -2708,7 +2918,7 @@ CONTAINS END IF END IF - ! always have the input param last to overwrite all the other ones + ! Always have the input param last to overwrite all the other ones IF (ASSOCIATED(inp_info%nonbonded14)) THEN DO k = 1, SIZE(inp_info%nonbonded14%pot) IF (iw > 0) WRITE (iw, *) " TESTING ", TRIM(name_atm_a), TRIM(name_atm_b), & @@ -2738,7 +2948,8 @@ CONTAINS END IF END DO END IF - ! at the very end we offer the possibility to overwrite the parameters for QM/MM + + ! At the very end we offer the possibility to overwrite the parameters for QM/MM ! nonbonded interactions IF (use_qmmm_ff) THEN match_names = 0 @@ -2747,9 +2958,14 @@ CONTAINS IF (match_names == 1) THEN IF (ASSOCIATED(qmmm_env%inp_info%nonbonded14)) THEN DO k = 1, SIZE(qmmm_env%inp_info%nonbonded14%pot) - IF (iw > 0) WRITE (iw, *) " TESTING ", TRIM(name_atm_a), TRIM(name_atm_b), & - " with ", TRIM(qmmm_env%inp_info%nonbonded14%pot(k)%pot%at1), & - TRIM(qmmm_env%inp_info%nonbonded14%pot(k)%pot%at2) + IF (debug_this_module) THEN + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(T2,A)") & + "FORCEFIELD| Testing "//TRIM(name_atm_a)//"-"//TRIM(name_atm_b)// & + " with "//TRIM(qmmm_env%inp_info%nonbonded14%pot(k)%pot%at1)//"-"// & + TRIM(qmmm_env%inp_info%nonbonded14%pot(k)%pot%at2) + END IF + END IF IF ((((name_atm_a) == (qmmm_env%inp_info%nonbonded14%pot(k)%pot%at1)) .AND. & ((name_atm_b) == (qmmm_env%inp_info%nonbonded14%pot(k)%pot%at2))) .OR. & (((name_atm_b) == (qmmm_env%inp_info%nonbonded14%pot(k)%pot%at1)) .AND. & @@ -2759,7 +2975,7 @@ CONTAINS IF (found) & CALL cp_warn(__LOCATION__, & "Multiple ONFO declaration: "//TRIM(name_atm_a)// & - " and "//TRIM(name_atm_b)//" ADDING QM/MM forcefield specifications! ") + " and "//TRIM(name_atm_b)//" Adding QM/MM forcefield specifications") potparm_nonbond14%pot(i, j)%pot => pot potparm_nonbond14%pot(j, i)%pot => pot ELSE @@ -2788,17 +3004,22 @@ CONTAINS pot%at1 = name_atm_a pot%at2 = name_atm_b END IF + ! If defined global RCUT let's use it IF (ff_type%rcut_nb > 0.0_dp) THEN pot%rcutsq = ff_type%rcut_nb*ff_type%rcut_nb END IF + ! Cutoff is defined always as the maximum between the FF and Ewald pot%rcutsq = MAX(pot%rcutsq, ewald_rcut*ewald_rcut) IF (only_qm) THEN CALL pair_potential_single_clean(pot) END IF + END DO ! atom kind j + END DO ! atom kind i + CALL timestop(handle2) END SUBROUTINE force_field_pack_nonbond14 @@ -2848,19 +3069,34 @@ CONTAINS TYPE(pair_potential_single_type), POINTER :: pot CALL timeset(routineN, handle2) + use_qmmm_ff = qmmm_env%use_qmmm_ff NULLIFY (pot) + + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(/,T2,A)") & + "FORCEFIELD| Checking for nonbonded terms" + END IF + CALL ewald_env_get(ewald_env, rcut=ewald_rcut) CALL pair_potential_pp_create(potparm_nonbond, SIZE(atomic_kind_set)) + DO i = 1, SIZE(atomic_kind_set) + atomic_kind => atomic_kind_set(i) + CALL get_atomic_kind(atomic_kind=atomic_kind, name=name_atm_a_local, & shell_active=is_a_shell) + DO j = i, SIZE(atomic_kind_set) + atomic_kind => atomic_kind_set(j) + CALL get_atomic_kind(atomic_kind=atomic_kind, name=name_atm_b_local, & shell_active=is_b_shell) + found = .FALSE. + name_atm_a = name_atm_a_local name_atm_b = name_atm_b_local only_qm = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b) @@ -2868,7 +3104,13 @@ CONTAINS CALL uppercase(name_atm_b) pot => potparm_nonbond%pot(i, j)%pot - ! loop over params from GROMOS + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(/,T2,A)") & + "FORCEFIELD| Checking for nonbonded terms between the atomic kinds "// & + TRIM(name_atm_a)//" and "//TRIM(name_atm_b) + END IF + + ! Loop over params from GROMOS IF (ASSOCIATED(gro_info%nonbond_a)) THEN ii = 0 jj = 0 @@ -2899,7 +3141,7 @@ CONTAINS END IF END IF - ! loop over params from CHARMM + ! Loop over params from CHARMM IF (ASSOCIATED(chm_info%nonbond_a)) THEN ii = 0 jj = 0 @@ -2931,7 +3173,7 @@ CONTAINS END IF END IF - ! loop over params from AMBER + ! Loop over params from AMBER IF (ASSOCIATED(amb_info%nonbond_a)) THEN ii = 0 jj = 0 @@ -2962,104 +3204,128 @@ CONTAINS END IF END IF - ! always have the input param last to overwrite all the other ones + ! Always have the input param last to overwrite all the other ones IF (ASSOCIATED(inp_info%nonbonded)) THEN DO k = 1, SIZE(inp_info%nonbonded%pot) IF ((TRIM(inp_info%nonbonded%pot(k)%pot%at1) == "*") .OR. & (TRIM(inp_info%nonbonded%pot(k)%pot%at2) == "*")) CYCLE - - IF (iw > 0) WRITE (iw, *) " TESTING ", TRIM(name_atm_a), TRIM(name_atm_b), & - " with ", TRIM(inp_info%nonbonded%pot(k)%pot%at1), & - TRIM(inp_info%nonbonded%pot(k)%pot%at2) + IF (debug_this_module) THEN + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(T2,A)") & + "FORCEFIELD| Testing "//TRIM(name_atm_a)//"-"//TRIM(name_atm_b)// & + " with "//TRIM(inp_info%nonbonded%pot(k)%pot%at1)//"-"// & + TRIM(inp_info%nonbonded%pot(k)%pot%at2) + END IF + END IF IF ((((name_atm_a) == (inp_info%nonbonded%pot(k)%pot%at1)) .AND. & ((name_atm_b) == (inp_info%nonbonded%pot(k)%pot%at2))) .OR. & (((name_atm_b) == (inp_info%nonbonded%pot(k)%pot%at1)) .AND. & ((name_atm_a) == (inp_info%nonbonded%pot(k)%pot%at2)))) THEN + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(T2,A)") & + "FORCEFIELD| Found nonbonded term "// & + TRIM(name_atm_a)//"-"//TRIM(name_atm_b) + END IF IF (ff_type%multiple_potential) THEN CALL pair_potential_single_add(inp_info%nonbonded%pot(k)%pot, pot) IF (found) & CALL cp_warn(__LOCATION__, & - "Multiple NONBONDED declaration: "//TRIM(name_atm_a)// & - " and "//TRIM(name_atm_b)//" ADDING! ") + "Multiple NONBONDED declarations for "//TRIM(name_atm_a)// & + "-"//TRIM(name_atm_b)//" -> ADDING") potparm_nonbond%pot(i, j)%pot => pot potparm_nonbond%pot(j, i)%pot => pot ELSE CALL pair_potential_single_copy(inp_info%nonbonded%pot(k)%pot, pot) IF (found) & CALL cp_warn(__LOCATION__, & - "Multiple NONBONDED declaration: "//TRIM(name_atm_a)// & - " and "//TRIM(name_atm_b)//" OVERWRITING! ") + "Multiple NONBONDED declarations for "//TRIM(name_atm_a)// & + "-"//TRIM(name_atm_b)//" -> OVERWRITING") END IF - IF (iw > 0) WRITE (iw, *) " FOUND ", TRIM(name_atm_a), " ", TRIM(name_atm_b) found = .TRUE. END IF END DO + ! Check for wildcards for one of the two types (if not associated yet) IF (.NOT. found) THEN DO k = 1, SIZE(inp_info%nonbonded%pot) IF ((TRIM(inp_info%nonbonded%pot(k)%pot%at1) == "*") .EQV. & (TRIM(inp_info%nonbonded%pot(k)%pot%at2) == "*")) CYCLE - - IF (iw > 0) WRITE (iw, *) " TESTING ", TRIM(name_atm_a), TRIM(name_atm_b), & - " with ", TRIM(inp_info%nonbonded%pot(k)%pot%at1), & - TRIM(inp_info%nonbonded%pot(k)%pot%at2) - + IF (debug_this_module) THEN + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(T2,A)") & + "FORCEFIELD| Testing "//TRIM(name_atm_a)//"-"//TRIM(name_atm_b)// & + " with "//TRIM(inp_info%nonbonded%pot(k)%pot%at1)//"-"// & + TRIM(inp_info%nonbonded%pot(k)%pot%at2) + END IF + END IF IF ((name_atm_a == inp_info%nonbonded%pot(k)%pot%at1) .OR. & (name_atm_b == inp_info%nonbonded%pot(k)%pot%at2) .OR. & (name_atm_b == inp_info%nonbonded%pot(k)%pot%at1) .OR. & (name_atm_a == inp_info%nonbonded%pot(k)%pot%at2)) THEN + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(T2,A)") & + "FORCEFIELD| Found one wildcard for "// & + TRIM(name_atm_a)//"-"//TRIM(name_atm_b) + END IF IF (ff_type%multiple_potential) THEN CALL pair_potential_single_add(inp_info%nonbonded%pot(k)%pot, pot) IF (found) & CALL cp_warn(__LOCATION__, & - "Multiple NONBONDED declaration: "//TRIM(name_atm_a)// & - " and "//TRIM(name_atm_b)//" ADDING! ") + "Multiple NONBONDED declarations "//TRIM(name_atm_a)// & + "-"//TRIM(name_atm_b)//" -> ADDING") potparm_nonbond%pot(i, j)%pot => pot potparm_nonbond%pot(j, i)%pot => pot ELSE CALL pair_potential_single_copy(inp_info%nonbonded%pot(k)%pot, pot) IF (found) & CALL cp_warn(__LOCATION__, & - "Multiple NONBONDED declaration: "//TRIM(name_atm_a)// & - " and "//TRIM(name_atm_b)//" OVERWRITING! ") + "Multiple NONBONDED declarations "//TRIM(name_atm_a)// & + "-"//TRIM(name_atm_b)//" -> OVERWRITING") END IF - IF (iw > 0) WRITE (iw, *) " FOUND (one WILDCARD)", TRIM(name_atm_a), " ", TRIM(name_atm_b) found = .TRUE. END IF END DO END IF + ! Check for wildcards for both types (if not associated yet) IF (.NOT. found) THEN DO k = 1, SIZE(inp_info%nonbonded%pot) IF ((TRIM(inp_info%nonbonded%pot(k)%pot%at1) /= "*") .OR. & (TRIM(inp_info%nonbonded%pot(k)%pot%at2) /= "*")) CYCLE - - IF (iw > 0) WRITE (iw, *) " TESTING ", TRIM(name_atm_a), TRIM(name_atm_b), & - " with ", TRIM(inp_info%nonbonded%pot(k)%pot%at1), & - TRIM(inp_info%nonbonded%pot(k)%pot%at2) - + IF (debug_this_module) THEN + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(T2,A)") & + "FORCEFIELD| Testing "//TRIM(name_atm_a)//"-"//TRIM(name_atm_b)// & + " with "//TRIM(inp_info%nonbonded%pot(k)%pot%at1)//"-"// & + TRIM(inp_info%nonbonded%pot(k)%pot%at2) + END IF + END IF + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(T2,A)") & + "FORCEFIELD| Found wildcards for both "// & + TRIM(name_atm_a)//" and "//TRIM(name_atm_b) + END IF IF (ff_type%multiple_potential) THEN CALL pair_potential_single_add(inp_info%nonbonded%pot(k)%pot, pot) IF (found) & CALL cp_warn(__LOCATION__, & - "Multiple NONBONDED declaration: "//TRIM(name_atm_a)// & - " and "//TRIM(name_atm_b)//" ADDING! ") + "Multiple NONBONDED declarations "//TRIM(name_atm_a)// & + " - "//TRIM(name_atm_b)//" -> ADDING") potparm_nonbond%pot(i, j)%pot => pot potparm_nonbond%pot(j, i)%pot => pot ELSE CALL pair_potential_single_copy(inp_info%nonbonded%pot(k)%pot, pot) IF (found) & CALL cp_warn(__LOCATION__, & - "Multiple NONBONDED declaration: "//TRIM(name_atm_a)// & - " and "//TRIM(name_atm_b)//" OVERWRITING! ") + "Multiple NONBONDED declarations "//TRIM(name_atm_a)// & + " - "//TRIM(name_atm_b)//" -> OVERWRITING") END IF - IF (iw > 0) WRITE (iw, *) " FOUND (both WILDCARDS)", TRIM(name_atm_a), " ", TRIM(name_atm_b) found = .TRUE. END DO END IF END IF - ! at the very end we offer the possibility to overwrite the parameters for QM/MM + ! At the very end we offer the possibility to overwrite the parameters for QM/MM ! nonbonded interactions IF (use_qmmm_ff) THEN match_names = 0 @@ -3068,35 +3334,44 @@ CONTAINS IF (match_names == 1) THEN IF (ASSOCIATED(qmmm_env%inp_info%nonbonded)) THEN DO k = 1, SIZE(qmmm_env%inp_info%nonbonded%pot) - IF (iw > 0) WRITE (iw, *) " TESTING ", TRIM(name_atm_a), TRIM(name_atm_b), & - " with ", TRIM(qmmm_env%inp_info%nonbonded%pot(k)%pot%at1), & - TRIM(qmmm_env%inp_info%nonbonded%pot(k)%pot%at2) + IF (debug_this_module) THEN + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(T2,A)") & + "FORCEFIELD| Testing "//TRIM(name_atm_a)//"-"//TRIM(name_atm_b)// & + " with "//TRIM(qmmm_env%inp_info%nonbonded%pot(k)%pot%at1), & + TRIM(qmmm_env%inp_info%nonbonded%pot(k)%pot%at2) + END IF + END IF IF ((((name_atm_a) == (qmmm_env%inp_info%nonbonded%pot(k)%pot%at1)) .AND. & ((name_atm_b) == (qmmm_env%inp_info%nonbonded%pot(k)%pot%at2))) .OR. & (((name_atm_b) == (qmmm_env%inp_info%nonbonded%pot(k)%pot%at1)) .AND. & ((name_atm_a) == (qmmm_env%inp_info%nonbonded%pot(k)%pot%at2)))) THEN + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(T2,A)") & + "FORCEFIELD| Found "//TRIM(name_atm_a)//"-"//TRIM(name_atm_b)//" (QM/MM)" + END IF IF (qmmm_env%multiple_potential) THEN CALL pair_potential_single_add(qmmm_env%inp_info%nonbonded%pot(k)%pot, pot) IF (found) & CALL cp_warn(__LOCATION__, & - "Multiple NONBONDED declaration: "//TRIM(name_atm_a)// & - " and "//TRIM(name_atm_b)//" ADDING QM/MM forcefield specifications! ") + "Multiple NONBONDED declarations for "//TRIM(name_atm_a)// & + " and "//TRIM(name_atm_b)//" -> ADDING QM/MM forcefield specifications") potparm_nonbond%pot(i, j)%pot => pot potparm_nonbond%pot(j, i)%pot => pot ELSE CALL pair_potential_single_copy(qmmm_env%inp_info%nonbonded%pot(k)%pot, pot) IF (found) & CALL cp_warn(__LOCATION__, & - "Multiple NONBONDED declaration: "//TRIM(name_atm_a)// & - " and "//TRIM(name_atm_b)//" OVERWRITING QM/MM forcefield specifications! ") + "Multiple NONBONDED declarations for "//TRIM(name_atm_a)// & + " and "//TRIM(name_atm_b)//" -> OVERWRITING QM/MM forcefield specifications") END IF - IF (iw > 0) WRITE (iw, *) " FOUND ", TRIM(name_atm_a), " ", TRIM(name_atm_b) found = .TRUE. END IF END DO END IF END IF END IF + IF (.NOT. found) THEN CALL store_FF_missing_par(atm1=TRIM(name_atm_a), & atm2=TRIM(name_atm_b), & @@ -3104,10 +3379,12 @@ CONTAINS fatal=fatal, & array=Ainfo) END IF + ! If defined global RCUT let's use it IF (ff_type%rcut_nb > 0.0_dp) THEN pot%rcutsq = ff_type%rcut_nb*ff_type%rcut_nb END IF + ! Cutoff is defined always as the maximum between the FF and Ewald pot%rcutsq = MAX(pot%rcutsq, ewald_rcut*ewald_rcut) ! Set the shell type @@ -3118,12 +3395,17 @@ CONTAINS ELSE pot%shell_type = nosh_nosh END IF + IF (only_qm) THEN CALL pair_potential_single_clean(pot) END IF + END DO ! jkind + END DO ! ikind + CALL timestop(handle2) + END SUBROUTINE force_field_pack_nonbond ! ************************************************************************************************** @@ -3154,20 +3436,27 @@ CONTAINS TYPE(spline_environment_type), POINTER :: spline_env CALL timeset(routineN, handle2) + + IF (iw2 > 0) THEN + WRITE (UNIT=iw2, FMT="(/,T2,A)") & + "FORCEFIELD| Splining nonbonded terms" + END IF + ! Figure out which nonbonded interactions happen to be identical, and ! prepare storage for these, avoiding duplicates. NULLIFY (spline_env) CALL get_nonbond_storage(spline_env, potparm, atomic_kind_set, & do_zbl, shift_cutoff=ff_type%shift_cutoff) - ! Effectively compute the spline data. + ! Effectively compute the spline data CALL spline_nonbond_control(spline_env, potparm, & atomic_kind_set, eps_spline=ff_type%eps_spline, & max_energy=ff_type%max_energy, rlow_nb=ff_type%rlow_nb, & - emax_spline=ff_type%emax_spline, npoints=ff_type%npoints, iw=iw2, iw2=iw3, iw3=iw4, & + emax_spline=ff_type%emax_spline, npoints=ff_type%npoints, & + iw=iw2, iw2=iw3, iw3=iw4, & do_zbl=do_zbl, shift_cutoff=ff_type%shift_cutoff, & nonbonded_type=nonbonded_type) ! Let the pointers on potparm point to the splines generated in - ! spline_nonbond_control. + ! spline_nonbond_control DO ikind = 1, SIZE(potparm%pot, 1) DO jkind = ikind, SIZE(potparm%pot, 2) n = spline_env%spltab(ikind, jkind) @@ -3180,6 +3469,12 @@ CONTAINS CALL spline_env_release(spline_env) DEALLOCATE (spline_env) NULLIFY (spline_env) + + IF (iw2 > 0) THEN + WRITE (UNIT=iw2, FMT="(/,T2,A)") & + "FORCEFIELD| Splining done" + END IF + CALL timestop(handle2) END SUBROUTINE force_field_pack_splines @@ -3190,15 +3485,16 @@ CONTAINS !> \param ff_type ... !> \param potparm_nonbond ... !> \param ewald_env ... +!> \param iw ... !> \author Toon.Verstraelen@gmail.com ! ************************************************************************************************** - SUBROUTINE force_field_pack_eicut(atomic_kind_set, ff_type, & - potparm_nonbond, ewald_env) + SUBROUTINE force_field_pack_eicut(atomic_kind_set, ff_type, potparm_nonbond, ewald_env, iw) TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(force_field_type), INTENT(IN) :: ff_type TYPE(pair_potential_pp_type), POINTER :: potparm_nonbond TYPE(ewald_environment_type), POINTER :: ewald_env + INTEGER, INTENT(IN) :: iw CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_eicut' @@ -3210,14 +3506,20 @@ CONTAINS CALL timeset(routineN, handle) + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(/,T2,A)") & + "FORCEFIELD| Computing the electrostatic interactions cutoffs" + END IF + tmp = 0.0_dp nkinds = SIZE(atomic_kind_set) - ! allocate the array with interaction cutoffs for the electrostatics, used + + ! Allocate the array with interaction cutoffs for the electrostatics, used ! to make the electrostatic interaction continuous at ewald_env%rcut ALLOCATE (interaction_cutoffs(3, nkinds, nkinds)) interaction_cutoffs = 0.0_dp - ! compute the interaction cutoff if SHIFT_CUTOFF is active + ! Compute the interaction cutoff if SHIFT_CUTOFF is active IF (ff_type%shift_cutoff) THEN CALL ewald_env_get(ewald_env, alpha=alpha, ewald_type=ewald_type, & rcut=rcut2_ewald) @@ -3260,6 +3562,7 @@ CONTAINS CALL ewald_env_set(ewald_env, interaction_cutoffs=interaction_cutoffs) CALL timestop(handle) + END SUBROUTINE force_field_pack_eicut ! ************************************************************************************************** @@ -3282,19 +3585,19 @@ CONTAINS CHARACTER(LEN=default_string_length) :: item - item = "( "//TRIM(name_atm_a) + item = "("//TRIM(name_atm_a) IF (PRESENT(name_atm_b)) THEN - item = TRIM(item)//" , "//TRIM(name_atm_b) + item = TRIM(item)//", "//TRIM(name_atm_b) END IF IF (PRESENT(name_atm_c)) THEN - item = TRIM(item)//" , "//TRIM(name_atm_c) + item = TRIM(item)//", "//TRIM(name_atm_c) END IF IF (PRESENT(name_atm_d)) THEN - item = TRIM(item)//" , "//TRIM(name_atm_d) + item = TRIM(item)//", "//TRIM(name_atm_d) END IF - item = TRIM(item)//" )" + item = TRIM(item)//")" IF (found) THEN - CPWARN("Multiple "//TRIM(tag_label)//" declarations: "//TRIM(item)//" overwriting!") + CPWARN("Found multiple "//TRIM(tag_label)//" terms for "//TRIM(item)//" -> OVERWRITING") END IF END SUBROUTINE issue_duplications @@ -3491,4 +3794,3 @@ CONTAINS END FUNCTION bsearch_leftmost_2d END MODULE force_fields_all - diff --git a/src/force_fields_util.F b/src/force_fields_util.F index 29b8af36ff..0315cc5935 100644 --- a/src/force_fields_util.F +++ b/src/force_fields_util.F @@ -164,86 +164,75 @@ CONTAINS gro_info => ff_type%gro_info amb_info => ff_type%amb_info !----------------------------------------------------------------------------- - !----------------------------------------------------------------------------- ! 1. Determine the number of unique bond kind and allocate bond_kind_set !----------------------------------------------------------------------------- CALL force_field_unique_bond(particle_set, molecule_kind_set, molecule_set, & - ff_type) - !----------------------------------------------------------------------------- + ff_type, iw) !----------------------------------------------------------------------------- ! 2. Determine the number of unique bend kind and allocate bend_kind_set !----------------------------------------------------------------------------- CALL force_field_unique_bend(particle_set, molecule_kind_set, molecule_set, & - ff_type) - !----------------------------------------------------------------------------- + ff_type, iw) !----------------------------------------------------------------------------- ! 3. Determine the number of unique Urey-Bradley kind and allocate ub_kind_set !----------------------------------------------------------------------------- - CALL force_field_unique_ub(particle_set, molecule_kind_set, molecule_set) - !----------------------------------------------------------------------------- + CALL force_field_unique_ub(particle_set, molecule_kind_set, molecule_set, iw) !----------------------------------------------------------------------------- ! 4. Determine the number of unique torsion kind and allocate torsion_kind_set !----------------------------------------------------------------------------- CALL force_field_unique_tors(particle_set, molecule_kind_set, molecule_set, & - ff_type) - !----------------------------------------------------------------------------- + ff_type, iw) !----------------------------------------------------------------------------- ! 5. Determine the number of unique impr kind and allocate impr_kind_set !----------------------------------------------------------------------------- CALL force_field_unique_impr(particle_set, molecule_kind_set, molecule_set, & - ff_type) - !----------------------------------------------------------------------------- + ff_type, iw) !----------------------------------------------------------------------------- ! 6. Determine the number of unique opbend kind and allocate opbend_kind_set !----------------------------------------------------------------------------- CALL force_field_unique_opbend(particle_set, molecule_kind_set, molecule_set, & - ff_type) - !----------------------------------------------------------------------------- + ff_type, iw) !----------------------------------------------------------------------------- ! 7. Bonds !----------------------------------------------------------------------------- CALL force_field_pack_bond(particle_set, molecule_kind_set, molecule_set, & - fatal, Ainfo, chm_info, inp_info, gro_info, amb_info) - !----------------------------------------------------------------------------- + fatal, Ainfo, chm_info, inp_info, gro_info, & + amb_info, iw) !----------------------------------------------------------------------------- ! 8. Bends !----------------------------------------------------------------------------- CALL force_field_pack_bend(particle_set, molecule_kind_set, molecule_set, & - fatal, Ainfo, chm_info, inp_info, gro_info, amb_info) + fatal, Ainfo, chm_info, inp_info, gro_info, & + amb_info, iw) ! Give information and abort if any bond or angle parameter is missing.. CALL release_FF_missing_par(fatal, ignore_fatal, AInfo, output_unit, iw) !----------------------------------------------------------------------------- - !----------------------------------------------------------------------------- ! 9. Urey-Bradley !----------------------------------------------------------------------------- CALL force_field_pack_ub(particle_set, molecule_kind_set, molecule_set, & Ainfo, chm_info, inp_info, iw) !----------------------------------------------------------------------------- - !----------------------------------------------------------------------------- ! 10. Torsions !----------------------------------------------------------------------------- CALL force_field_pack_tors(particle_set, molecule_kind_set, molecule_set, & Ainfo, chm_info, inp_info, gro_info, amb_info, iw) !----------------------------------------------------------------------------- - !----------------------------------------------------------------------------- ! 11. Impropers !----------------------------------------------------------------------------- CALL force_field_pack_impr(particle_set, molecule_kind_set, molecule_set, & - Ainfo, chm_info, inp_info, gro_info) - !----------------------------------------------------------------------------- + Ainfo, chm_info, inp_info, gro_info, iw) !----------------------------------------------------------------------------- ! 12. Out of plane bends !----------------------------------------------------------------------------- CALL force_field_pack_opbend(particle_set, molecule_kind_set, molecule_set, & - Ainfo, inp_info) + Ainfo, inp_info, iw) ! Give information only if any Urey-Bradley, Torsion, improper or opbend is missing - ! continue calculation.. + ! continue calculation CALL release_FF_missing_par(fatal, ignore_fatal, AInfo, output_unit, iw) charges_section => section_vals_get_subs_vals(mm_section, "FORCEFIELD%CHARGES") CALL section_vals_get(charges_section, explicit=explicit) IF (.NOT. explicit) THEN - !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- ! 13.a Set up atomic_kind_set()%fist_potentail%[qeff] and shell ! potential parameters @@ -253,7 +242,6 @@ CONTAINS ! Give information only if charge is missing and abort.. CALL release_FF_missing_par(fatal, ignore_fatal, AInfo, output_unit, iw) ELSE - !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- ! 13.b Setup a global array of classical charges - this avoids the packing and ! allows the usage of different charges for same atomic types @@ -261,21 +249,15 @@ CONTAINS CALL force_field_pack_charges(charges, charges_section, particle_set, my_qmmm, & qmmm_env, inp_info, iw4) END IF - - !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- ! 14. Set up the radius of the electrostatic multipole in Fist !----------------------------------------------------------------------------- CALL force_field_pack_radius(atomic_kind_set, iw, subsys_section) - - !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- ! 15. Set up the polarizable FF parameters !----------------------------------------------------------------------------- CALL force_field_pack_pol(atomic_kind_set, iw, inp_info) CALL force_field_pack_damp(atomic_kind_set, iw, inp_info) - - !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- ! 16. Set up Shell potential !----------------------------------------------------------------------------- @@ -283,7 +265,6 @@ CONTAINS molecule_kind_set, molecule_set, root_section, subsys_section, & shell_particle_set, core_particle_set, cell, iw, inp_info) IF (ff_type%do_nonbonded) THEN - !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- ! 17. Set up potparm_nonbond14 !----------------------------------------------------------------------------- @@ -297,7 +278,6 @@ CONTAINS CALL force_field_pack_splines(atomic_kind_set, ff_type, iw2, iw3, iw4, & potparm_nonbond14, do_zbl, nonbonded_type="NONBONDED14") !----------------------------------------------------------------------------- - !----------------------------------------------------------------------------- ! 18. Set up potparm_nonbond !----------------------------------------------------------------------------- ! Move the data from the info structures to potparm_nonbond @@ -312,7 +292,6 @@ CONTAINS potparm_nonbond, do_zbl, nonbonded_type="NONBONDED") END IF !----------------------------------------------------------------------------- - !----------------------------------------------------------------------------- ! 19. Create nonbond environment !----------------------------------------------------------------------------- CALL ewald_env_get(ewald_env, rcut=ewald_rcut) @@ -326,12 +305,13 @@ CONTAINS CALL fist_nonbond_env_set(fist_nonbond_env, charges=charges) ! Compute the electrostatic interaction cutoffs. CALL force_field_pack_eicut(atomic_kind_set, ff_type, potparm_nonbond, & - ewald_env) + ewald_env, iw) CALL cp_print_key_finished_output(iw4, logger, mm_section, "PRINT%PROGRAM_RUN_INFO") CALL cp_print_key_finished_output(iw3, logger, mm_section, "PRINT%FF_INFO/SPLINE_DATA") CALL cp_print_key_finished_output(iw2, logger, mm_section, "PRINT%FF_INFO/SPLINE_INFO") CALL cp_print_key_finished_output(iw, logger, mm_section, "PRINT%FF_INFO") + CALL timestop(handle) END SUBROUTINE force_field_pack @@ -420,6 +400,7 @@ CONTAINS TYPE(shell_kind_type), POINTER :: shell CALL timeset(routineN, handle) + NULLIFY (logger) logger => cp_get_default_logger() iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_INFO", & @@ -429,8 +410,6 @@ CONTAINS qeff_mol = 0.0_dp qeff_sum = 0.0_dp mass_sum = 0.0_dp - - !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- ! 1. Sum of [qeff,mass] for each molecule_kind !----------------------------------------------------------------------------- @@ -461,8 +440,6 @@ CONTAINS CALL set_molecule_kind(molecule_kind=molecule_kind, charge=qeff_mol, mass=mass_mol) IF (iw > 0) WRITE (iw, *) " Mol Kind ", TRIM(molname), " charge = ", qeff_mol END DO - - !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- ! 2. Sum of [qeff,mass] for particle_set !----------------------------------------------------------------------------- @@ -484,7 +461,9 @@ CONTAINS CALL cp_print_key_finished_output(iw, logger, mm_section, & "PRINT%FF_INFO") + CALL timestop(handle) + END SUBROUTINE force_field_qeff_output ! ************************************************************************************************** @@ -529,7 +508,6 @@ CONTAINS iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_INFO", & extension=".mmLog") !----------------------------------------------------------------------------- - !----------------------------------------------------------------------------- ! 1. Lets Tag the unwanted bonds due to the use of distance constraint !----------------------------------------------------------------------------- DO ikind = 1, SIZE(molecule_kind_set) @@ -557,7 +535,6 @@ CONTAINS END IF END DO !----------------------------------------------------------------------------- - !----------------------------------------------------------------------------- ! 2. Lets Tag the unwanted bends due to the use of distance constraint !----------------------------------------------------------------------------- DO ikind = 1, SIZE(molecule_kind_set) @@ -591,7 +568,6 @@ CONTAINS END IF END DO !----------------------------------------------------------------------------- - !----------------------------------------------------------------------------- ! 3. Lets Tag the unwanted bonds due to the use of 3x3 !----------------------------------------------------------------------------- DO ikind = 1, SIZE(molecule_kind_set) @@ -617,7 +593,6 @@ CONTAINS END DO END DO !----------------------------------------------------------------------------- - !----------------------------------------------------------------------------- ! 4. Lets Tag the unwanted bends due to the use of 3x3 !----------------------------------------------------------------------------- DO ikind = 1, SIZE(molecule_kind_set) @@ -645,7 +620,6 @@ CONTAINS END DO END DO !----------------------------------------------------------------------------- - !----------------------------------------------------------------------------- ! 5. Lets Tag the unwanted bonds due to the use of 4x6 !----------------------------------------------------------------------------- DO ikind = 1, SIZE(molecule_kind_set) @@ -672,7 +646,6 @@ CONTAINS END DO END DO !----------------------------------------------------------------------------- - !----------------------------------------------------------------------------- ! 6. Lets Tag the unwanted bends due to the use of 4x6 !----------------------------------------------------------------------------- DO ikind = 1, SIZE(molecule_kind_set) @@ -700,7 +673,6 @@ CONTAINS END DO END DO !----------------------------------------------------------------------------- - !----------------------------------------------------------------------------- ! 7. Count the number of UNSET bond kinds there are !----------------------------------------------------------------------------- DO ikind = 1, SIZE(molecule_kind_set) @@ -782,7 +754,6 @@ CONTAINS END IF END DO !----------------------------------------------------------------------------- - !----------------------------------------------------------------------------- ! 8. Count the number of UNSET bend kinds there are !----------------------------------------------------------------------------- DO ikind = 1, SIZE(molecule_kind_set) @@ -864,8 +835,6 @@ CONTAINS DEALLOCATE (bad1) END IF END DO - - !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- ! 9. Count the number of UNSET Urey-Bradley kinds there are !----------------------------------------------------------------------------- @@ -948,8 +917,6 @@ CONTAINS DEALLOCATE (bad1) END IF END DO - - !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- ! 10. Count the number of UNSET torsion kinds there are !----------------------------------------------------------------------------- @@ -1045,8 +1012,6 @@ CONTAINS DEALLOCATE (bad1) END IF END DO - - !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- ! 11. Count the number of UNSET improper kinds there are !----------------------------------------------------------------------------- @@ -1128,8 +1093,6 @@ CONTAINS DEALLOCATE (bad1) END IF END DO - - !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- ! 11. Count the number of UNSET opbends kinds there are !----------------------------------------------------------------------------- @@ -1208,15 +1171,15 @@ CONTAINS DEALLOCATE (bad1) END IF END DO - - !--------------------------------------------------------------------------- !--------------------------------------------------------------------------- ! 12. Count the number of UNSET NONBOND14 kinds there are !- NEED TO REMOVE EXTRAS HERE - IKUO !--------------------------------------------------------------------------- CALL cp_print_key_finished_output(iw, logger, mm_section, & "PRINT%FF_INFO") + CALL timestop(handle) + END SUBROUTINE clean_intra_force_kind ! ************************************************************************************************** diff --git a/src/nnp_acsf.F b/src/nnp_acsf.F index d3a88f4fc1..ad0ec8cf0d 100644 --- a/src/nnp_acsf.F +++ b/src/nnp_acsf.F @@ -71,7 +71,7 @@ CONTAINS REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: forcetmp REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: force3tmp REAL(KIND=dp), DIMENSION(3) :: rvect1, rvect2, rvect3 - TYPE(nnp_neighbor_type) :: neighbor + TYPE(nnp_neighbor_type) :: neighbor = nnp_neighbor_type() CALL timeset(routineN, handle) diff --git a/src/pair_potential.F b/src/pair_potential.F index 204ec19133..56be9156b9 100644 --- a/src/pair_potential.F +++ b/src/pair_potential.F @@ -132,9 +132,9 @@ CONTAINS !> \par History !> Teo 2006.05 : Improved speed and accuracy. Linear scaling of the setup ! ************************************************************************************************** - SUBROUTINE spline_nonbond_control(spline_env, potparm, atomic_kind_set, & - eps_spline, max_energy, rlow_nb, emax_spline, npoints, iw, iw2, iw3, do_zbl, & - shift_cutoff, nonbonded_type) + SUBROUTINE spline_nonbond_control(spline_env, potparm, atomic_kind_set, eps_spline, & + max_energy, rlow_nb, emax_spline, npoints, iw, iw2, iw3, & + do_zbl, shift_cutoff, nonbonded_type) TYPE(spline_environment_type), POINTER :: spline_env TYPE(pair_potential_pp_type), POINTER :: potparm @@ -153,11 +153,12 @@ CONTAINS REAL(KIND=dp) :: energy_cutoff, hicut, hicut0, locut TYPE(pair_potential_single_type), POINTER :: pot + CALL timeset(routineN, handle) + n = 0 ncount = 0 - ntype = SIZE(atomic_kind_set) - CALL timeset(routineN, handle) + IF (iw3 > 0) THEN WRITE (iw3, "(/,T2,A,I0,A,I0,A)") & "SPLINE_INFO| Generating ", (ntype*(ntype + 1))/2, " splines for "// & @@ -255,13 +256,12 @@ CONTAINS CALL finalizef() IF (iw > 0) THEN - WRITE (iw, '(/,T2,A,I0)') & + WRITE (UNIT=iw, FMT='(/,T2,A,I0)') & "SPLINE_INFO| Number of pair potential splines allocated: ", MAXVAL(spline_env%spltab) END IF IF (iw3 > 0) THEN - WRITE (iw3, '(/,T2,A,I0,/,T2,A)') & - "SPLINE_INFO| Number of unique splines computed: ", MAXVAL(spline_env%spltab), & - "SPLINE_INFO| Done" + WRITE (UNIT=iw3, FMT='(/,T2,A,I0)') & + "SPLINE_INFO| Number of unique splines computed: ", MAXVAL(spline_env%spltab) END IF CALL timestop(handle) diff --git a/src/qs_rho0_methods.F b/src/qs_rho0_methods.F index 451aa72968..c64d9000db 100644 --- a/src/qs_rho0_methods.F +++ b/src/qs_rho0_methods.F @@ -492,7 +492,9 @@ CONTAINS IF (ASSOCIATED(cneo_potential)) THEN IF (PRESENT(zcore)) THEN - IF (zcore == 0.0_dp) CPABORT("Electronic TDDFT with CNEO quantum nuclei is not implemented.") + IF (zcore == 0.0_dp) THEN + CPABORT("Electronic TDDFT with CNEO quantum nuclei is not implemented.") + END IF END IF CPASSERT(paw_atom) NULLIFY (nuc_basis, nuc_soft_basis) diff --git a/src/subsys/molecule_kind_types.F b/src/subsys/molecule_kind_types.F index 2eb2596a29..5ffaae2b10 100644 --- a/src/subsys/molecule_kind_types.F +++ b/src/subsys/molecule_kind_types.F @@ -11,7 +11,7 @@ !> \par History !> Teodoro Laino [tlaino] 12.2008 - Preparing for VIRTUAL SITE constraints !> (patch by Marcel Baer) -!> \author MK (22.08.2003) +!> \author Matthias Krack (22.08.2003) ! ************************************************************************************************** MODULE molecule_kind_types USE atomic_kind_types, ONLY: atomic_kind_type,& @@ -50,159 +50,153 @@ MODULE molecule_kind_types CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'molecule_kind_types' -! *** Define the derived structure types *** + ! Define the derived structure types -! ************************************************************************************************** TYPE atom_type TYPE(atomic_kind_type), POINTER :: atomic_kind => NULL() INTEGER :: id_name = 0 END TYPE atom_type -! ************************************************************************************************** TYPE shell_type INTEGER :: a = 0 - CHARACTER(LEN=default_string_length) :: name = "" + CHARACTER(LEN=default_string_length) :: name = "" TYPE(shell_kind_type), POINTER :: shell_kind => NULL() END TYPE shell_type -! ************************************************************************************************** TYPE bond_type INTEGER :: a = 0, b = 0 INTEGER :: id_type = do_ff_undef, itype = 0 TYPE(bond_kind_type), POINTER :: bond_kind => NULL() END TYPE bond_type -! ************************************************************************************************** TYPE bend_type INTEGER :: a = 0, b = 0, c = 0 INTEGER :: id_type = do_ff_undef, itype = 0 TYPE(bend_kind_type), POINTER :: bend_kind => NULL() END TYPE bend_type -! ************************************************************************************************** TYPE ub_type INTEGER :: a = 0, b = 0, c = 0 INTEGER :: id_type = do_ff_undef, itype = 0 TYPE(ub_kind_type), POINTER :: ub_kind => NULL() END TYPE ub_type -! ************************************************************************************************** TYPE torsion_type INTEGER :: a = 0, b = 0, c = 0, d = 0 INTEGER :: id_type = do_ff_undef, itype = 0 TYPE(torsion_kind_type), POINTER :: torsion_kind => NULL() END TYPE torsion_type -! ************************************************************************************************** TYPE impr_type INTEGER :: a = 0, b = 0, c = 0, d = 0 INTEGER :: id_type = do_ff_undef, itype = 0 TYPE(impr_kind_type), POINTER :: impr_kind => NULL() END TYPE impr_type -! ************************************************************************************************** TYPE opbend_type INTEGER :: a = 0, b = 0, c = 0, d = 0 INTEGER :: id_type = do_ff_undef, itype = 0 TYPE(opbend_kind_type), POINTER :: opbend_kind => NULL() END TYPE opbend_type -! ************************************************************************************************** TYPE restraint_type - LOGICAL :: active = .FALSE. + LOGICAL :: active = .FALSE. REAL(KIND=dp) :: k0 = 0.0_dp END TYPE restraint_type -! ************************************************************************************************** + ! Constraint types TYPE colvar_constraint_type INTEGER :: type_id = no_colvar_id INTEGER :: inp_seq_num = 0 LOGICAL :: use_points = .FALSE. - REAL(KIND=dp) :: expected_value = 0.0_dp - REAL(KIND=dp) :: expected_value_growth_speed = 0.0_dp + REAL(KIND=dp) :: expected_value = 0.0_dp + REAL(KIND=dp) :: expected_value_growth_speed = 0.0_dp INTEGER, POINTER, DIMENSION(:) :: i_atoms => NULL() TYPE(restraint_type) :: restraint = restraint_type() END TYPE colvar_constraint_type -! ************************************************************************************************** TYPE g3x3_constraint_type INTEGER :: a = 0, b = 0, c = 0 - REAL(KIND=dp) :: dab = 0.0_dp, dac = 0.0_dp, dbc = 0.0_dp + REAL(KIND=dp) :: dab = 0.0_dp, dac = 0.0_dp, dbc = 0.0_dp TYPE(restraint_type) :: restraint = restraint_type() END TYPE g3x3_constraint_type -! ************************************************************************************************** TYPE g4x6_constraint_type INTEGER :: a = 0, b = 0, c = 0, d = 0 - REAL(KIND=dp) :: dab = 0.0_dp, dac = 0.0_dp, dbc = 0.0_dp, dad = 0.0_dp, dbd = 0.0_dp, dcd = 0.0_dp + REAL(KIND=dp) :: dab = 0.0_dp, dac = 0.0_dp, dbc = 0.0_dp, & + dad = 0.0_dp, dbd = 0.0_dp, dcd = 0.0_dp TYPE(restraint_type) :: restraint = restraint_type() END TYPE g4x6_constraint_type -! ************************************************************************************************** TYPE vsite_constraint_type INTEGER :: a = 0, b = 0, c = 0, d = 0 - REAL(KIND=dp) :: wbc = 0.0_dp, wdc = 0.0_dp + REAL(KIND=dp) :: wbc = 0.0_dp, wdc = 0.0_dp TYPE(restraint_type) :: restraint = restraint_type() END TYPE vsite_constraint_type -! ************************************************************************************************** TYPE fixd_constraint_type TYPE(restraint_type) :: restraint = restraint_type() INTEGER :: fixd = 0, itype = 0 REAL(KIND=dp), DIMENSION(3) :: coord = 0.0_dp END TYPE fixd_constraint_type -! ************************************************************************************************** TYPE local_fixd_constraint_type INTEGER :: ifixd_index = 0, ikind = 0 END TYPE local_fixd_constraint_type -! ************************************************************************************************** + ! Molecule kind type TYPE molecule_kind_type - TYPE(atom_type), DIMENSION(:), POINTER :: atom_list => NULL() - TYPE(bond_kind_type), DIMENSION(:), POINTER :: bond_kind_set => NULL() - TYPE(bond_type), DIMENSION(:), POINTER :: bond_list => NULL() - TYPE(bend_kind_type), DIMENSION(:), POINTER :: bend_kind_set => NULL() - TYPE(bend_type), DIMENSION(:), POINTER :: bend_list => NULL() - TYPE(ub_kind_type), DIMENSION(:), POINTER :: ub_kind_set => NULL() - TYPE(ub_type), DIMENSION(:), POINTER :: ub_list => NULL() - TYPE(torsion_kind_type), DIMENSION(:), POINTER :: torsion_kind_set => NULL() - TYPE(torsion_type), DIMENSION(:), POINTER :: torsion_list => NULL() - TYPE(impr_kind_type), DIMENSION(:), POINTER :: impr_kind_set => NULL() - TYPE(impr_type), DIMENSION(:), POINTER :: impr_list => NULL() - TYPE(opbend_kind_type), DIMENSION(:), POINTER :: opbend_kind_set => NULL() - TYPE(opbend_type), DIMENSION(:), POINTER :: opbend_list => NULL() - TYPE(colvar_constraint_type), DIMENSION(:), POINTER :: colv_list => NULL() - TYPE(g3x3_constraint_type), DIMENSION(:), POINTER :: g3x3_list => NULL() - TYPE(g4x6_constraint_type), DIMENSION(:), POINTER :: g4x6_list => NULL() - TYPE(vsite_constraint_type), DIMENSION(:), POINTER :: vsite_list => NULL() - TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list => NULL() - TYPE(shell_type), DIMENSION(:), POINTER :: shell_list => NULL() - CHARACTER(LEN=default_string_length) :: name = "" - REAL(KIND=dp) :: charge = 0.0_dp, & - mass = 0.0_dp - INTEGER :: kind_number = 0, & - natom = 0, & - nbond = 0, & - nbend = 0, & - nimpr = 0, & - nopbend = 0, & - ntorsion = 0, & - nub = 0, & - ng3x3 = 0, ng3x3_restraint = 0, & - ng4x6 = 0, ng4x6_restraint = 0, & - nvsite = 0, nvsite_restraint = 0, & - nfixd = 0, nfixd_restraint = 0, & - nmolecule = 0, nshell = 0 - TYPE(colvar_counters) :: ncolv = colvar_counters() - INTEGER :: nsgf = 0, nelectron = 0, & - nelectron_alpha = 0, & - nelectron_beta = 0 - INTEGER, DIMENSION(:), POINTER :: molecule_list => NULL() - LOGICAL :: molname_generated = .FALSE. + TYPE(atom_type), DIMENSION(:), POINTER :: atom_list => NULL() + TYPE(bond_kind_type), DIMENSION(:), POINTER :: bond_kind_set => NULL() + TYPE(bond_type), DIMENSION(:), POINTER :: bond_list => NULL() + TYPE(bend_kind_type), DIMENSION(:), POINTER :: bend_kind_set => NULL() + TYPE(bend_type), DIMENSION(:), POINTER :: bend_list => NULL() + TYPE(ub_kind_type), DIMENSION(:), POINTER :: ub_kind_set => NULL() + TYPE(ub_type), DIMENSION(:), POINTER :: ub_list => NULL() + TYPE(torsion_kind_type), DIMENSION(:), POINTER :: torsion_kind_set => NULL() + TYPE(torsion_type), DIMENSION(:), POINTER :: torsion_list => NULL() + TYPE(impr_kind_type), DIMENSION(:), POINTER :: impr_kind_set => NULL() + TYPE(impr_type), DIMENSION(:), POINTER :: impr_list => NULL() + TYPE(opbend_kind_type), DIMENSION(:), POINTER :: opbend_kind_set => NULL() + TYPE(opbend_type), DIMENSION(:), POINTER :: opbend_list => NULL() + TYPE(colvar_constraint_type), DIMENSION(:), & + POINTER :: colv_list => NULL() + TYPE(g3x3_constraint_type), DIMENSION(:), POINTER :: g3x3_list => NULL() + TYPE(g4x6_constraint_type), DIMENSION(:), POINTER :: g4x6_list => NULL() + TYPE(vsite_constraint_type), DIMENSION(:), POINTER :: vsite_list => NULL() + TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list => NULL() + TYPE(shell_type), DIMENSION(:), POINTER :: shell_list => NULL() + CHARACTER(LEN=default_string_length) :: name = "" + REAL(KIND=dp) :: charge = 0.0_dp, & + mass = 0.0_dp + INTEGER :: kind_number = 0, & + natom = 0, & + nbond = 0, & + nbend = 0, & + nimpr = 0, & + nopbend = 0, & + ntorsion = 0, & + nub = 0, & + ng3x3 = 0, & + ng3x3_restraint = 0, & + ng4x6 = 0, & + ng4x6_restraint = 0, & + nvsite = 0, & + nvsite_restraint = 0, & + nfixd = 0, & + nfixd_restraint = 0, & + nmolecule = 0, & + nshell = 0 + TYPE(colvar_counters) :: ncolv = colvar_counters() + INTEGER :: nsgf = 0, & + nelectron = 0, & + nelectron_alpha = 0, & + nelectron_beta = 0 + INTEGER, DIMENSION(:), POINTER :: molecule_list => NULL() + LOGICAL :: molname_generated = .FALSE. END TYPE molecule_kind_type - ! *** Public subroutines *** + ! Public subroutines PUBLIC :: allocate_molecule_kind_set, & deallocate_molecule_kind_set, & get_molecule_kind, & @@ -211,7 +205,7 @@ MODULE molecule_kind_types write_molecule_kind_set, & setup_colvar_counters - ! *** Public data types *** + ! Public data types PUBLIC :: atom_type, & bend_type, & bond_type, & @@ -316,7 +310,7 @@ CONTAINS !> \param molecule_kind_set ... !> \param nmolecule_kind ... !> \date 22.08.2003 -!> \author MK +!> \author Matthias Krack !> \version 1.0 ! ************************************************************************************************** SUBROUTINE allocate_molecule_kind_set(molecule_kind_set, nmolecule_kind) @@ -343,7 +337,7 @@ CONTAINS !> \brief Deallocate a molecule kind set. !> \param molecule_kind_set ... !> \date 22.08.2003 -!> \author MK +!> \author Matthias Krack !> \version 1.0 ! ************************************************************************************************** SUBROUTINE deallocate_molecule_kind_set(molecule_kind_set) @@ -493,7 +487,7 @@ CONTAINS !> \param torsion_kind_set ... !> \param molname_generated ... !> \date 27.08.2003 -!> \author MK +!> \author Matthias Krack !> \version 1.0 ! ************************************************************************************************** SUBROUTINE get_molecule_kind(molecule_kind, atom_list, bond_list, bend_list, & @@ -648,7 +642,7 @@ CONTAINS !> \param nmolecule ... !> \param nrestraints ... !> \date 27.08.2003 -!> \author MK +!> \author Matthias Krack !> \version 1.0 ! ************************************************************************************************** SUBROUTINE get_molecule_kind_set(molecule_kind_set, maxatom, natom, & @@ -761,7 +755,7 @@ CONTAINS !> \param nsgf ... !> \param molname_generated ... !> \date 27.08.2003 -!> \author MK +!> \author Matthias Krack !> \version 1.0 ! ************************************************************************************************** SUBROUTINE set_molecule_kind(molecule_kind, name, mass, charge, kind_number, & @@ -889,7 +883,7 @@ CONTAINS !> \param molecule_kind ... !> \param output_unit ... !> \date 24.09.2003 -!> \author MK +!> \author Matthias Krack !> \version 1.0 ! ************************************************************************************************** SUBROUTINE write_molecule_kind(molecule_kind, output_unit) @@ -959,7 +953,7 @@ CONTAINS !> \param molecule_kind_set ... !> \param subsys_section ... !> \date 24.09.2003 -!> \author MK +!> \author Matthias Krack !> \version 1.0 ! ************************************************************************************************** SUBROUTINE write_molecule_kind_set(molecule_kind_set, subsys_section) diff --git a/src/subsys/molecule_types.F b/src/subsys/molecule_types.F index f6c48b8677..61776084f3 100644 --- a/src/subsys/molecule_types.F +++ b/src/subsys/molecule_types.F @@ -11,7 +11,7 @@ !> JGH (22.05.2004) add last_atom information !> Teodoro Laino [tlaino] 12.2008 - Preparing for VIRTUAL SITE constraints !> (patch by Marcel Baer) -!> \author MK (29.08.2003) +!> \author Matthias Krack (29.08.2003) ! ************************************************************************************************** MODULE molecule_types @@ -32,92 +32,141 @@ MODULE molecule_types PRIVATE -! *** Global parameters (in this module) *** + ! Global parameters (in this module) CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'molecule_types' -! *** Data types *** -! ************************************************************************************************** + ! Molecular constraint types TYPE local_colvar_constraint_type - LOGICAL :: init = .FALSE. - TYPE(colvar_type), POINTER :: colvar => NULL() - TYPE(colvar_type), POINTER :: colvar_old => NULL() - REAL(KIND=dp) :: lambda = 0.0_dp, sigma = 0.0_dp + TYPE(colvar_type), POINTER :: colvar => NULL(), & + colvar_old => NULL() + REAL(KIND=dp) :: lambda = 0.0_dp, & + sigma = 0.0_dp + LOGICAL :: init = .FALSE. END TYPE local_colvar_constraint_type -! ************************************************************************************************** TYPE local_g3x3_constraint_type - LOGICAL :: init = .FALSE. - REAL(KIND=dp) :: scale = 0.0_dp, scale_old = 0.0_dp, & - imass1 = 0.0_dp, imass2 = 0.0_dp, imass3 = 0.0_dp - REAL(KIND=dp), DIMENSION(3) :: fa = 0.0_dp, fb = 0.0_dp, fc = 0.0_dp, & - f_roll1 = 0.0_dp, f_roll2 = 0.0_dp, f_roll3 = 0.0_dp, & - ra_old = 0.0_dp, rb_old = 0.0_dp, rc_old = 0.0_dp, & - va = 0.0_dp, vb = 0.0_dp, vc = 0.0_dp, & - lambda = 0.0_dp, del_lambda = 0.0_dp, lambda_old = 0.0_dp, & - r0_12 = 0.0_dp, r0_13 = 0.0_dp, r0_23 = 0.0_dp - REAL(KIND=dp), DIMENSION(3, 3) :: amat = 0.0_dp + LOGICAL :: init = .FALSE. + REAL(KIND=dp) :: scale = 0.0_dp, & + imass1 = 0.0_dp, & + imass2 = 0.0_dp, & + imass3 = 0.0_dp, & + scale_old = 0.0_dp + REAL(KIND=dp), DIMENSION(3) :: fa = 0.0_dp, & + fb = 0.0_dp, & + fc = 0.0_dp, & + f_roll1 = 0.0_dp, & + f_roll2 = 0.0_dp, & + f_roll3 = 0.0_dp, & + ra_old = 0.0_dp, & + rb_old = 0.0_dp, & + rc_old = 0.0_dp, & + r0_12 = 0.0_dp, & + r0_13 = 0.0_dp, & + r0_23 = 0.0_dp, & + va = 0.0_dp, & + vb = 0.0_dp, & + vc = 0.0_dp, & + del_lambda = 0.0_dp, & + lambda = 0.0_dp, & + lambda_old = 0.0_dp + REAL(KIND=dp), DIMENSION(3, 3) :: amat = 0.0_dp END TYPE local_g3x3_constraint_type -! ************************************************************************************************** TYPE local_g4x6_constraint_type - LOGICAL :: init = .FALSE. - REAL(KIND=dp) :: scale = 0.0_dp, scale_old = 0.0_dp, imass1 = 0.0_dp, & - imass2 = 0.0_dp, imass3 = 0.0_dp, imass4 = 0.0_dp - REAL(KIND=dp), DIMENSION(3) :: fa = 0.0_dp, fb = 0.0_dp, fc = 0.0_dp, fd = 0.0_dp, fe = 0.0_dp, ff = 0.0_dp, & - f_roll1 = 0.0_dp, f_roll2 = 0.0_dp, f_roll3 = 0.0_dp, & - f_roll4 = 0.0_dp, f_roll5 = 0.0_dp, f_roll6 = 0.0_dp, & - ra_old = 0.0_dp, rb_old = 0.0_dp, rc_old = 0.0_dp, & - rd_old = 0.0_dp, re_old = 0.0_dp, rf_old = 0.0_dp, & - va = 0.0_dp, vb = 0.0_dp, vc = 0.0_dp, vd = 0.0_dp, ve = 0.0_dp, vf = 0.0_dp, & - r0_12 = 0.0_dp, r0_13 = 0.0_dp, r0_14 = 0.0_dp, & - r0_23 = 0.0_dp, r0_24 = 0.0_dp, r0_34 = 0.0_dp - REAL(KIND=dp), DIMENSION(6) :: lambda = 0.0_dp, del_lambda = 0.0_dp, lambda_old = 0.0_dp - REAL(KIND=dp), DIMENSION(6, 6) :: amat = 0.0_dp + LOGICAL :: init = .FALSE. + REAL(KIND=dp) :: scale = 0.0_dp, & + scale_old = 0.0_dp, & + imass1 = 0.0_dp, & + imass2 = 0.0_dp, & + imass3 = 0.0_dp, & + imass4 = 0.0_dp + REAL(KIND=dp), DIMENSION(3) :: fa = 0.0_dp, & + fb = 0.0_dp, & + fc = 0.0_dp, & + fd = 0.0_dp, & + fe = 0.0_dp, & + ff = 0.0_dp, & + f_roll1 = 0.0_dp, & + f_roll2 = 0.0_dp, & + f_roll3 = 0.0_dp, & + f_roll4 = 0.0_dp, & + f_roll5 = 0.0_dp, & + f_roll6 = 0.0_dp, & + ra_old = 0.0_dp, & + rb_old = 0.0_dp, & + rc_old = 0.0_dp, & + rd_old = 0.0_dp, & + re_old = 0.0_dp, & + rf_old = 0.0_dp, & + va = 0.0_dp, & + vb = 0.0_dp, & + vc = 0.0_dp, & + vd = 0.0_dp, & + ve = 0.0_dp, & + vf = 0.0_dp, & + r0_12 = 0.0_dp, & + r0_13 = 0.0_dp, & + r0_14 = 0.0_dp, & + r0_23 = 0.0_dp, & + r0_24 = 0.0_dp, & + r0_34 = 0.0_dp + REAL(KIND=dp), DIMENSION(6) :: del_lambda = 0.0_dp, & + lambda = 0.0_dp, & + lambda_old = 0.0_dp + REAL(KIND=dp), DIMENSION(6, 6) :: amat = 0.0_dp END TYPE local_g4x6_constraint_type -! ************************************************************************************************** TYPE local_states_type - INTEGER, DIMENSION(:), POINTER :: states => NULL() ! indices of Kohn-Sham states for molecule - INTEGER :: nstates = 0 ! Kohn-Sham states for molecule + INTEGER :: nstates = 0 ! Kohn-Sham states for molecule + INTEGER, DIMENSION(:), POINTER :: states => NULL() ! indices of Kohn-Sham states for molecule END TYPE local_states_type -! ************************************************************************************************** TYPE local_constraint_type - TYPE(local_colvar_constraint_type), DIMENSION(:), POINTER :: lcolv => NULL() - TYPE(local_g3x3_constraint_type), DIMENSION(:), POINTER :: lg3x3 => NULL() - TYPE(local_g4x6_constraint_type), DIMENSION(:), POINTER :: lg4x6 => NULL() + TYPE(local_colvar_constraint_type), & + DIMENSION(:), POINTER :: lcolv => NULL() + TYPE(local_g3x3_constraint_type), DIMENSION(:), & + POINTER :: lg3x3 => NULL() + TYPE(local_g4x6_constraint_type), DIMENSION(:), & + POINTER :: lg4x6 => NULL() END TYPE local_constraint_type -! ************************************************************************************************** TYPE global_constraint_type - TYPE(colvar_counters) :: ncolv = colvar_counters() - INTEGER :: ntot = 0, nrestraint = 0 - INTEGER :: ng3x3 = 0, ng3x3_restraint = 0 - INTEGER :: ng4x6 = 0, ng4x6_restraint = 0 - INTEGER :: nvsite = 0, nvsite_restraint = 0 - TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list => NULL() - TYPE(colvar_constraint_type), DIMENSION(:), POINTER :: colv_list => NULL() - TYPE(g3x3_constraint_type), DIMENSION(:), POINTER :: g3x3_list => NULL() - TYPE(g4x6_constraint_type), DIMENSION(:), POINTER :: g4x6_list => NULL() - TYPE(vsite_constraint_type), DIMENSION(:), POINTER :: vsite_list => NULL() - TYPE(local_colvar_constraint_type), DIMENSION(:), POINTER :: lcolv => NULL() - TYPE(local_g3x3_constraint_type), DIMENSION(:), POINTER :: lg3x3 => NULL() - TYPE(local_g4x6_constraint_type), DIMENSION(:), POINTER :: lg4x6 => NULL() + TYPE(colvar_counters) :: ncolv = colvar_counters() + INTEGER :: ntot = 0, & + nrestraint = 0, & + ng3x3 = 0, & + ng3x3_restraint = 0, & + ng4x6 = 0, & + ng4x6_restraint = 0, & + nvsite = 0, & + nvsite_restraint = 0 + TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list => NULL() + TYPE(colvar_constraint_type), DIMENSION(:), & + POINTER :: colv_list => NULL() + TYPE(g3x3_constraint_type), DIMENSION(:), POINTER :: g3x3_list => NULL() + TYPE(g4x6_constraint_type), DIMENSION(:), POINTER :: g4x6_list => NULL() + TYPE(vsite_constraint_type), DIMENSION(:), POINTER :: vsite_list => NULL() + TYPE(local_colvar_constraint_type), & + DIMENSION(:), POINTER :: lcolv => NULL() + TYPE(local_g3x3_constraint_type), DIMENSION(:), & + POINTER :: lg3x3 => NULL() + TYPE(local_g4x6_constraint_type), DIMENSION(:), & + POINTER :: lg4x6 => NULL() END TYPE global_constraint_type -! ************************************************************************************************** + ! Molecule type TYPE molecule_type - TYPE(molecule_kind_type), POINTER :: molecule_kind => NULL() ! pointer to molecule kind information - TYPE(local_states_type), DIMENSION(:), POINTER :: lmi => NULL() ! local (spin)-states information - TYPE(local_constraint_type), POINTER :: lci => NULL() ! local molecule constraint info - INTEGER :: first_atom = 0 ! global index of first atom in molecule - INTEGER :: last_atom = 0 ! global index of last atom in molecule - INTEGER :: first_shell = 0 ! global index of first shell atom in molecule - INTEGER :: last_shell = 0 ! global index of last shell atom in molecule + TYPE(molecule_kind_type), POINTER :: molecule_kind => NULL() ! pointer to molecule kind information + TYPE(local_states_type), DIMENSION(:), POINTER :: lmi => NULL() ! local (spin)-states information + TYPE(local_constraint_type), POINTER :: lci => NULL() ! local molecule constraint info + INTEGER :: first_atom = 0 ! global index of first atom in molecule + INTEGER :: last_atom = 0 ! global index of last atom in molecule + INTEGER :: first_shell = 0 ! global index of first shell atom in molecule + INTEGER :: last_shell = 0 ! global index of last shell atom in molecule END TYPE molecule_type -! *** Public data types *** + ! Public data types PUBLIC :: local_colvar_constraint_type, & local_g3x3_constraint_type, & @@ -127,7 +176,7 @@ MODULE molecule_types global_constraint_type, & molecule_type -! *** Public subroutines *** + ! Public subroutines PUBLIC :: deallocate_global_constraint, & allocate_molecule_set, & @@ -195,7 +244,7 @@ CONTAINS !> \param molecule_set ... !> \param nmolecule ... !> \date 29.08.2003 -!> \author MK +!> \author Matthias Krack !> \version 1.0 ! ************************************************************************************************** SUBROUTINE allocate_molecule_set(molecule_set, nmolecule) @@ -212,7 +261,7 @@ CONTAINS !> \brief Deallocate a molecule set. !> \param molecule_set ... !> \date 29.08.2003 -!> \author MK +!> \author Matthias Krack !> \version 1.0 ! ************************************************************************************************** SUBROUTINE deallocate_molecule_set(molecule_set) @@ -269,7 +318,7 @@ CONTAINS !> \param first_shell ... !> \param last_shell ... !> \date 29.08.2003 -!> \author MK +!> \author Matthias Krack !> \version 1.0 ! ************************************************************************************************** SUBROUTINE get_molecule(molecule, molecule_kind, lmi, lci, lg3x3, lg4x6, lcolv, & @@ -330,7 +379,7 @@ CONTAINS !> \param lg3x3 ... !> \param lg4x6 ... !> \date 29.08.2003 -!> \author MK +!> \author Matthias Krack !> \version 1.0 ! ************************************************************************************************** SUBROUTINE set_molecule(molecule, molecule_kind, lmi, lci, lcolv, lg3x3, lg4x6) @@ -379,7 +428,7 @@ CONTAINS !> \param first_atom ... !> \param last_atom ... !> \date 29.08.2003 -!> \author MK +!> \author Matthias Krack !> \version 1.0 ! ************************************************************************************************** SUBROUTINE set_molecule_set(molecule_set, first_atom, last_atom) diff --git a/src/topology.F b/src/topology.F index 5a6082bf43..9a271a6f57 100644 --- a/src/topology.F +++ b/src/topology.F @@ -97,7 +97,7 @@ MODULE topology PRIVATE -! *** Public parameters *** + ! Public parameters PUBLIC :: topology_control, & connectivity_control @@ -538,8 +538,6 @@ CONTAINS NULLIFY (global_section) global_section => section_vals_get_subs_vals(root_section, "GLOBAL") CALL section_vals_val_get(global_section, "SAVE_MEM", l_val=save_mem) - - !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- ! 1. If reading in from external file, make sure its there first !----------------------------------------------------------------------------- @@ -558,7 +556,6 @@ CONTAINS END IF END IF !----------------------------------------------------------------------------- - !----------------------------------------------------------------------------- ! 2. Read in the coordinates from the corresponding locations !----------------------------------------------------------------------------- CALL timeset(routineN//"_READ_COORDINATE", handle2) diff --git a/src/topology_generate_util.F b/src/topology_generate_util.F index 75aeb3ec57..bc5d2c5826 100644 --- a/src/topology_generate_util.F +++ b/src/topology_generate_util.F @@ -99,34 +99,35 @@ CONTAINS CHARACTER(LEN=default_string_length), PARAMETER :: basename = "MOL" CHARACTER(LEN=default_string_length) :: molname - INTEGER :: i, N, nmol + INTEGER :: i, id_undef, n, nmol LOGICAL :: check TYPE(array1_list_type), ALLOCATABLE, DIMENSION(:) :: atom_bond_list -! convert a simple list of bonds to a list of bonds per atom -! (each bond is present in the forward and backward direction + ! convert a simple list of bonds to a list of bonds per atom + ! (each bond is present in the forward and backward direction ALLOCATE (atom_bond_list(natom)) - DO I = 1, natom - ALLOCATE (atom_bond_list(I)%array1(0)) + DO i = 1, natom + ALLOCATE (atom_bond_list(i)%array1(0)) END DO - N = 0 - IF (ASSOCIATED(conn_info%bond_a)) N = SIZE(conn_info%bond_a) - nbond_prev + n = 0 + IF (ASSOCIATED(conn_info%bond_a)) n = SIZE(conn_info%bond_a) - nbond_prev CALL reorder_structure(atom_bond_list, conn_info%bond_a(nbond_prev + 1:) - natom_prev, & - conn_info%bond_b(nbond_prev + 1:) - natom_prev, N) + conn_info%bond_b(nbond_prev + 1:) - natom_prev, n) nmol = 0 - check = ALL(id_molname == str2id(s2s("__UNDEF__"))) .OR. ALL(id_molname /= str2id(s2s("__UNDEF__"))) + id_undef = str2id(s2s("__UNDEF__")) + check = ALL(id_molname == id_undef) .OR. ALL(id_molname /= id_undef) CPASSERT(check) DO i = 1, natom - IF (id2str(id_molname(i)) == "__UNDEF__") THEN + IF (id_molname(i) == id_undef) THEN molname = TRIM(basename)//ADJUSTL(cp_to_string(nmol)) CALL generate_molname_low(i, atom_bond_list, molname, id_molname) nmol = nmol + 1 END IF END DO - DO I = 1, natom - DEALLOCATE (atom_bond_list(I)%array1) + DO i = 1, natom + DEALLOCATE (atom_bond_list(i)%array1) END DO DEALLOCATE (atom_bond_list) @@ -243,7 +244,8 @@ CONTAINS ELSE IF (topology%conn_type == do_conn_user) THEN ! User-defined connectivity: 5th column of COORD section or molecule ! or residue name in the case of PDB files - IF (atom_info%id_molname(iatom) == atom_info%id_molname(iatom - 1)) THEN + IF ((atom_info%id_molname(iatom) == atom_info%id_molname(iatom - 1)) .AND. & + (.NOT. MODULO(iatom, topology%natom_muc) == 1)) THEN atom_info%map_mol_typ(iatom) = atom_info%map_mol_typ(iatom - 1) IF (atom_info%id_resname(iatom) == atom_info%id_resname(iatom - 1)) THEN atom_info%map_mol_res(iatom) = atom_info%map_mol_res(iatom - 1) diff --git a/src/topology_multiple_unit_cell.F b/src/topology_multiple_unit_cell.F index b73d55d9ff..097a914103 100644 --- a/src/topology_multiple_unit_cell.F +++ b/src/topology_multiple_unit_cell.F @@ -29,13 +29,13 @@ MODULE topology_multiple_unit_cell PRIVATE -! *** Public parameters *** + ! Public parameters PUBLIC :: topology_muc CONTAINS ! ************************************************************************************************** -!> \brief Handles the multiple_unit_cell for the atomic coordinates.. +!> \brief Handles the multiple_unit_cell for the atomic coordinates. !> \param topology ... !> \param subsys_section ... !> \author Teodoro Laino [tlaino] - 05.2009 @@ -55,21 +55,30 @@ CONTAINS TYPE(section_vals_type), POINTER :: work_section CALL timeset(routineN, handle) + NULLIFY (multiple_unit_cell, iwork, cell) + + ! Store original number of atoms for the molecule generation in any case + topology%natom_muc = topology%natoms + CALL section_vals_val_get(subsys_section, "TOPOLOGY%MULTIPLE_UNIT_CELL", & i_vals=multiple_unit_cell) + ! Fail is one of the value is set to zero.. IF (ANY(multiple_unit_cell <= 0)) & CALL cp_abort(__LOCATION__, "SUBSYS%TOPOLOGY%MULTIPLE_UNIT_CELL accepts "// & - "only integer values larger than 0! A value of 0 or negative is meaningless!") + "only integer values greater than zero.") + IF (ANY(multiple_unit_cell /= 1)) THEN - ! Check that the setup between CELL and TOPOLOGY is the same.. + + ! Check that the setup between CELL and TOPOLOGY is the same CALL section_vals_val_get(subsys_section, "CELL%MULTIPLE_UNIT_CELL", & i_vals=iwork) IF (ANY(iwork /= multiple_unit_cell)) & - CALL cp_abort(__LOCATION__, "SUBSYS%TOPOLOGY%MULTIPLE_UNIT_CELL and "// & - "SUBSYS%CELL%MULTIPLE_UNIT_CELL have been "// & - "setup to two different values!! Correct this error!") + CALL cp_abort(__LOCATION__, "The input parameters for "// & + "SUBSYS%TOPOLOGY%MULTIPLE_UNIT_CELL and "// & + "SUBSYS%CELL%MULTIPLE_UNIT_CELL have to agree.") + cell => topology%cell_muc natoms = topology%natoms*PRODUCT(multiple_unit_cell) @@ -80,8 +89,8 @@ CONTAINS CALL section_vals_val_get(work_section, '_DEFAULT_KEYWORD_', n_rep_val=nrep) check = nrep == natoms IF (.NOT. check) & - CALL cp_abort(__LOCATION__, & - "Number of available entries in VELOCITY section is not compatible with the number of atoms!") + CALL cp_abort(__LOCATION__, "The number of available entries in the "// & + "VELOCITY section is not compatible with the number of atoms.") END IF CALL reallocate(topology%atom_info%id_molname, 1, natoms) @@ -94,6 +103,7 @@ CONTAINS CALL reallocate(topology%atom_info%occup, 1, natoms) CALL reallocate(topology%atom_info%beta, 1, natoms) CALL reallocate(topology%atom_info%id_element, 1, natoms) + ind = 0 DO k = 1, multiple_unit_cell(3) trsl_k = cell%hmat(:, 3)*REAL(k - 1, KIND=dp) @@ -104,8 +114,7 @@ CONTAINS trsl = trsl_i + trsl_j + trsl_k ind = ind + 1 IF (ind == 1) CYCLE - - ! loop over atoms + ! Loop over all atoms n = (ind - 1)*topology%natoms DO m = 1, topology%natoms topology%atom_info%id_atmname(n + m) = topology%atom_info%id_atmname(m) @@ -122,6 +131,7 @@ CONTAINS END DO END DO END DO + ! Store the new total number of atoms topology%natoms = natoms ! Deallocate the coordinate section (will be rebuilt later with the whole atomic set) @@ -137,7 +147,9 @@ CONTAINS CALL section_vals_val_set(work_section, "SCALED", l_val=scale) END IF END IF + CALL timestop(handle) + END SUBROUTINE topology_muc END MODULE topology_multiple_unit_cell diff --git a/src/topology_types.F b/src/topology_types.F index 4e00a94463..f47b436460 100644 --- a/src/topology_types.F +++ b/src/topology_types.F @@ -33,20 +33,20 @@ MODULE topology_types ! ************************************************************************************************** TYPE atom_info_type - INTEGER, DIMENSION(:), POINTER :: id_molname => NULL() - INTEGER, DIMENSION(:), POINTER :: id_resname => NULL() - INTEGER, DIMENSION(:), POINTER :: id_atmname => NULL() - INTEGER, DIMENSION(:), POINTER :: id_atom_names => NULL() - INTEGER, DIMENSION(:), POINTER :: id_element => NULL() - INTEGER, POINTER :: resid(:) => NULL() - REAL(KIND=dp), DIMENSION(:, :), POINTER :: r => NULL() - INTEGER, POINTER :: map_mol_typ(:) => NULL() - INTEGER, POINTER :: map_mol_num(:) => NULL() - INTEGER, POINTER :: map_mol_res(:) => NULL() - REAL(KIND=dp), POINTER :: atm_charge(:) => NULL() - REAL(KIND=dp), POINTER :: atm_mass(:) => NULL() - REAL(KIND=dp), POINTER :: occup(:) => NULL() - REAL(KIND=dp), POINTER :: beta(:) => NULL() + INTEGER, DIMENSION(:), POINTER :: id_molname => NULL(), & + id_resname => NULL(), & + id_atmname => NULL(), & + id_atom_names => NULL(), & + id_element => NULL() + INTEGER, POINTER :: resid(:) => NULL() + REAL(KIND=dp), DIMENSION(:, :), POINTER :: r => NULL() + INTEGER, POINTER :: map_mol_typ(:) => NULL() + INTEGER, POINTER :: map_mol_num(:) => NULL() + INTEGER, POINTER :: map_mol_res(:) => NULL() + REAL(KIND=dp), POINTER :: atm_charge(:) => NULL() + REAL(KIND=dp), POINTER :: atm_mass(:) => NULL() + REAL(KIND=dp), POINTER :: occup(:) => NULL() + REAL(KIND=dp), POINTER :: beta(:) => NULL() END TYPE atom_info_type ! ************************************************************************************************** @@ -144,14 +144,14 @@ MODULE topology_types TYPE(atom_info_type), POINTER :: atom_info => NULL() TYPE(connectivity_info_type), POINTER :: conn_info => NULL() TYPE(constraint_info_type), POINTER :: cons_info => NULL() - TYPE(cell_type), POINTER :: cell => NULL(), cell_ref => NULL(), cell_muc => NULL() + TYPE(cell_type), POINTER :: cell => NULL(), cell_ref => NULL(), cell_muc => NULL() INTEGER :: conn_type = -1 INTEGER :: coord_type = -1 INTEGER :: exclude_vdw = -1 INTEGER :: exclude_ei = -1 INTEGER :: bondparm_type = -1 !TRY TO REMOVE THIS FIVE VARIABLE IN THE FUTURE - INTEGER :: natoms = -1, natom_type = -1 + INTEGER :: natoms = -1, natom_type = -1, natom_muc = -1 INTEGER :: nmol = -1, nmol_type = -1, nmol_conn = -1 !TRY TO REMOVE THIS FIVE VARIABLE IN THE FUTURE LOGICAL :: aa_element = .FALSE. @@ -162,8 +162,8 @@ MODULE topology_types LOGICAL :: molecules_check = .FALSE. LOGICAL :: coordinate = .FALSE. LOGICAL :: use_g96_velocity = .FALSE. - CHARACTER(LEN=default_path_length) :: coord_file_name = "" - CHARACTER(LEN=default_path_length) :: conn_file_name = "" + CHARACTER(LEN=default_path_length) :: coord_file_name = "" + CHARACTER(LEN=default_path_length) :: conn_file_name = "" LOGICAL :: const_atom = .FALSE. LOGICAL :: const_hydr = .FALSE. LOGICAL :: const_colv = .FALSE. @@ -217,6 +217,7 @@ CONTAINS !----------------------------------------------------------------------------- NULLIFY (topology%cell, topology%cell_ref, topology%cell_muc) topology%natoms = 0 + topology%natom_muc = 0 topology%natom_type = 0 topology%nmol = 0 topology%nmol_type = 0 diff --git a/src/topology_util.F b/src/topology_util.F index 8a7d161e5f..ceaae74b1f 100644 --- a/src/topology_util.F +++ b/src/topology_util.F @@ -46,12 +46,10 @@ MODULE topology_util CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_util' -! ************************************************************************************************** TYPE array1_list_type INTEGER, POINTER, DIMENSION(:) :: array1 => NULL() END TYPE array1_list_type -! ************************************************************************************************** TYPE array2_list_type INTEGER, POINTER, DIMENSION(:) :: array1 => NULL() INTEGER, POINTER, DIMENSION(:) :: array2 => NULL() @@ -1009,20 +1007,24 @@ CONTAINS TYPE(connectivity_info_type), POINTER :: conn_info TYPE(cp_logger_type), POINTER :: logger + CALL timeset(routineN, handle) + NULLIFY (logger) logger => cp_get_default_logger() iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/UTIL_INFO", & extension=".subsysLog") - CALL timeset(routineN, handle) atom_info => topology%atom_info conn_info => topology%conn_info natom = topology%natoms - IF (iw > 0) WRITE (iw, '(A)') "Start of Molecule_Check", & - " Checking consistency between the generated molecules" + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(/,T2,A)") & + "FORCEFIELD| Checking consistency of the generated molecules" + END IF ALLOCATE (atom_bond_list(natom)) + DO I = 1, natom ALLOCATE (atom_bond_list(I)%array1(0)) END DO @@ -1040,17 +1042,18 @@ CONTAINS icheck_num = (atom_info%map_mol_num(iatom) == mol_num) icheck_typ = (atom_info%map_mol_typ(iatom) == mol_typ) IF ((icheck_typ .AND. (.NOT. icheck_num)) .OR. (.NOT. icheck_typ)) THEN - !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- ! 1. Check each molecule have the same number of atoms !----------------------------------------------------------------------------- IF (counter /= loc_counter) THEN CALL cp_abort(__LOCATION__, & - "different number of atoms for same molecule kind"// & - " molecule type = "//cp_to_string(mol_typ)// & - " molecule number= "//cp_to_string(mol_num)// & - " expected number of atoms="//cp_to_string(counter)//" found="// & - cp_to_string(loc_counter)) + "Different number of atoms for the same molecule kind found "// & + "(molecule type = "//TRIM(ADJUSTL(cp_to_string(mol_typ)))// & + ", molecule number = "//TRIM(ADJUSTL(cp_to_string(mol_num)))// & + ", expected number of atoms = "// & + TRIM(ADJUSTL(cp_to_string(counter)))// & + ", number of atoms found = "// & + TRIM(ADJUSTL(cp_to_string(loc_counter)))//").") END IF END IF IF (.NOT. icheck_typ) THEN @@ -1063,42 +1066,42 @@ CONTAINS IF (icheck_num) THEN IF (icheck_typ) loc_counter = loc_counter + 1 !----------------------------------------------------------------------------- - !----------------------------------------------------------------------------- ! 2. Check that each molecule has the same atom sequences !----------------------------------------------------------------------------- IF (atom_info%id_atmname(iatom) /= & atom_info%id_atmname(first + loc_counter - 1)) THEN CALL cp_abort(__LOCATION__, & - "different atom name for same molecule kind"// & - " atom number = "//cp_to_string(iatom)// & - " molecule type = "//cp_to_string(mol_typ)// & - " molecule number= "//cp_to_string(mol_num)// & - " expected atom name="//TRIM(id2str(atom_info%id_atmname(first + loc_counter - 1)))// & - " found="//TRIM(id2str(atom_info%id_atmname(iatom)))) + "Different atom names for the same molecule kind found "// & + "(atom number = "//TRIM(ADJUSTL(cp_to_string(iatom)))// & + ", molecule type = "//TRIM(ADJUSTL(cp_to_string(mol_typ)))// & + ", molecule number = "//TRIM(ADJUSTL(cp_to_string(mol_num)))// & + ", expected atom name = "// & + TRIM(ADJUSTL(id2str(atom_info%id_atmname(first + loc_counter - 1))))// & + ", atom name found = "// & + TRIM(ADJUSTL(id2str(atom_info%id_atmname(iatom))))//").") END IF !----------------------------------------------------------------------------- - !----------------------------------------------------------------------------- ! 3. Check that each molecule have the same bond sequences !----------------------------------------------------------------------------- IF (SIZE(atom_bond_list(iatom)%array1) /= SIZE(atom_bond_list(first + loc_counter - 1)%array1)) THEN CALL cp_abort(__LOCATION__, & - "different number of bonds for same molecule kind"// & - " molecule type = "//cp_to_string(mol_typ)// & - " molecule number= "//cp_to_string(mol_num)// & - " expected bonds="// & - cp_to_string(SIZE(atom_bond_list(first + loc_counter - 1)%array1))//" - "// & - cp_to_string(SIZE(atom_bond_list(iatom)%array1))// & - " NOT FOUND! Check the connectivity of your system.") + "Different number of bonds for the same molecule kind found "// & + "(molecule type = "//TRIM(ADJUSTL(cp_to_string(mol_typ)))// & + ", molecule number = "//TRIM(ADJUSTL(cp_to_string(mol_num)))// & + ", expected number of bonds = "// & + TRIM(ADJUSTL(cp_to_string(SIZE(atom_bond_list(first + loc_counter - 1)%array1))))// & + ", number of bonds found = "// & + TRIM(ADJUSTL(cp_to_string(SIZE(atom_bond_list(iatom)%array1))))// & + "). Check the connectivity of your system.") END IF - DO k = 1, SIZE(atom_bond_list(iatom)%array1) IF (ALL(atom_bond_list(first + loc_counter - 1)%array1 - first /= & atom_bond_list(iatom)%array1(k) - first_loc)) THEN CALL cp_abort(__LOCATION__, & - "different sequence of bonds for same molecule kind"// & - " molecule type = "//cp_to_string(mol_typ)// & - " molecule number= "//cp_to_string(mol_num)// & - " NOT FOUND! Check the connectivity of your system.") + "Different sequence of bonds for the same molecule kind found "// & + "(molecule type = "//TRIM(ADJUSTL(cp_to_string(mol_typ)))// & + ", molecule number = "//TRIM(ADJUSTL(cp_to_string(mol_num)))// & + "). Check the connectivity of your system.") END IF END DO ELSE @@ -1108,16 +1111,22 @@ CONTAINS END IF IF (mol_num == 1 .AND. icheck_typ) counter = counter + 1 END DO - IF (iw > 0) WRITE (iw, '(A)') "End of Molecule_Check" + + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(T2,A)") & + "FORCEFIELD| Consistency check for molecules completed" + END IF DO I = 1, natom DEALLOCATE (atom_bond_list(I)%array1) END DO DEALLOCATE (atom_bond_list) - CALL timestop(handle) + CALL cp_print_key_finished_output(iw, logger, subsys_section, & "PRINT%TOPOLOGY_INFO/UTIL_INFO") + CALL timestop(handle) + END SUBROUTINE topology_molecules_check ! **************************************************************************************************