comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-28 22:25:32 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Fix molecule indexing when using MULTIPLE_UNIT_CELL

Browse source 
Revise forcefield output
This commit is contained in:
Matthias Krack 2025-09-12 19:32:01 +02:00
parent 8bdb955e8a
commit 0446f17dce
12 changed files with 787 additions and 456 deletions
Download patch file Download diff file Expand all files Collapse all files

632
src/force_fields_all.F
View file

File diff suppressed because it is too large Load diff

77
src/force_fields_util.F
View file

@ -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
! **************************************************************************************************

2
src/nnp_acsf.F
View file

@ -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)

18
src/pair_potential.F
View file

@ -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)

4
src/qs_rho0_methods.F
View file

@ -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)

144
src/subsys/molecule_kind_types.F
View file

@ -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)

187
src/subsys/molecule_types.F
View file

@ -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)

5
src/topology.F
View file

@ -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)

28
src/topology_generate_util.F
View file

@ -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)

34
src/topology_multiple_unit_cell.F
View file

@ -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

37
src/topology_types.F
View file

@ -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

75
src/topology_util.F
View file

@ -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
! **************************************************************************************************