From c83fa4e215c08dfa6c3e5f1582f0d0c82b64389b Mon Sep 17 00:00:00 2001 From: cschran Date: Tue, 20 Aug 2019 17:05:46 +0200 Subject: [PATCH] Added centroid constraints for PINT runs --- src/motion/pint_methods.F | 524 +++++++++++++++++- src/motion/pint_types.F | 20 +- src/start/input_cp2k_motion.F | 6 + tests/Pimd/regtest-1/TEST_FILES | 2 + .../regtest-1/h2o_pint_constraints_exact.inp | 67 +++ .../regtest-1/h2o_pint_constraints_nose.inp | 69 +++ 6 files changed, 657 insertions(+), 31 deletions(-) create mode 100644 tests/Pimd/regtest-1/h2o_pint_constraints_exact.inp create mode 100644 tests/Pimd/regtest-1/h2o_pint_constraints_nose.inp diff --git a/src/motion/pint_methods.F b/src/motion/pint_methods.F index 485f9127d4..1aba64d6f7 100644 --- a/src/motion/pint_methods.F +++ b/src/motion/pint_methods.F @@ -15,11 +15,18 @@ ! ************************************************************************************************** MODULE pint_methods - USE atomic_kind_types, ONLY: get_atomic_kind + USE atomic_kind_list_types, ONLY: atomic_kind_list_type + USE atomic_kind_types, ONLY: atomic_kind_type,& + get_atomic_kind USE bibliography, ONLY: Brieuc2016,& Ceriotti2010,& Ceriotti2012,& cite_reference + USE cell_types, ONLY: cell_type + USE constraint, ONLY: rattle_control,& + shake_control,& + shake_update_targets + USE constraint_util, ONLY: getold USE cp_external_control, ONLY: external_control USE cp_log_handling, ONLY: cp_get_default_logger,& cp_logger_get_default_io_unit,& @@ -27,7 +34,9 @@ MODULE pint_methods cp_to_string USE cp_output_handling, ONLY: cp_add_iter_level,& cp_iterate,& + cp_p_file,& cp_print_key_finished_output,& + cp_print_key_should_output,& cp_print_key_unit_nr,& cp_rm_iter_level USE cp_para_types, ONLY: cp_para_env_type @@ -35,6 +44,7 @@ MODULE pint_methods cp_subsys_type USE cp_units, ONLY: cp_unit_from_cp2k,& cp_unit_to_cp2k + USE distribution_1d_types, ONLY: distribution_1d_type USE f77_interface, ONLY: f_env_add_defaults,& f_env_rm_defaults,& f_env_type @@ -71,7 +81,13 @@ MODULE pint_methods m_walltime USE mathconstants, ONLY: twopi USE mathlib, ONLY: gcd - USE message_passing, ONLY: mp_bcast + USE message_passing, ONLY: mp_bcast,& + mp_comm_self + USE molecule_kind_list_types, ONLY: molecule_kind_list_type + USE molecule_kind_types, ONLY: molecule_kind_type + USE molecule_list_types, ONLY: molecule_list_type + USE molecule_types, ONLY: global_constraint_type,& + molecule_type USE parallel_rng_types, ONLY: GAUSSIAN,& create_rng_stream,& delete_rng_stream,& @@ -79,6 +95,7 @@ MODULE pint_methods next_rng_seed,& rng_stream_type USE particle_list_types, ONLY: particle_list_type + USE particle_types, ONLY: particle_type USE pint_gle, ONLY: pint_calc_gle_energy,& pint_gle_init,& pint_gle_step @@ -123,6 +140,7 @@ MODULE pint_methods rep_env_create USE replica_types, ONLY: rep_env_release,& replica_env_type + USE simpar_types, ONLY: create_simpar_type #include "../base/base_uses.f90" IMPLICIT NONE @@ -146,6 +164,7 @@ CONTAINS !> Fixed some bugs [hforbert] !> Added normal mode transformation [hforbert] !> 10.2015 Added RPMD propagator and harmonic integrator [Felix Uhl] +!> 10.2018 Added centroid constraints [cschran+rperez] !> \author fawzi !> \note Might return an unassociated pointer in parallel on the processors !> that are not needed. @@ -161,23 +180,22 @@ CONTAINS CHARACTER(len=2*default_string_length) :: msg CHARACTER(len=default_path_length) :: output_file_name, project_name - INTEGER :: handle, iat, ibead, idim, idir, ierr, & - itmp, nrep, prep, stat + INTEGER :: handle, iat, ibead, icont, idim, idir, & + ierr, ig, itmp, nrep, prep, stat LOGICAL :: explicit, ltmp REAL(kind=dp) :: dt, mass, omega REAL(kind=dp), DIMENSION(3, 2) :: seed TYPE(cp_subsys_type), POINTER :: subsys TYPE(f_env_type), POINTER :: f_env + TYPE(global_constraint_type), POINTER :: gci TYPE(particle_list_type), POINTER :: particles TYPE(replica_env_type), POINTER :: rep_env - TYPE(section_vals_type), POINTER :: gle_section, nose_section, & - piglet_section, pile_section, & - pint_section, qtb_section, & - transform_section + TYPE(section_vals_type), POINTER :: constraint_section, gle_section, nose_section, & + piglet_section, pile_section, pint_section, qtb_section, transform_section CALL timeset(routineN, handle) - NULLIFY (f_env, subsys, particles, nose_section, gle_section) + NULLIFY (f_env, subsys, particles, nose_section, gle_section, gci) CPASSERT(.NOT. ASSOCIATED(pint_env)) CPASSERT(ASSOCIATED(input)) @@ -401,7 +419,7 @@ CONTAINS CALL f_env_add_defaults(f_env_id=pint_env%replicas%f_env_id, & f_env=f_env) CALL force_env_get(force_env=f_env%force_env, subsys=subsys) - CALL cp_subsys_get(subsys, particles=particles) + CALL cp_subsys_get(subsys, particles=particles, gci=gci) !TODO length of pint_env%mass is redundant idim = 0 @@ -595,6 +613,79 @@ CONTAINS pint_env%first_propagated_mode = 1 END IF + ! Constraint information: + NULLIFY (pint_env%simpar) + constraint_section => section_vals_get_subs_vals(pint_env%input, & + "MOTION%CONSTRAINT") + CALL section_vals_get(constraint_section, explicit=explicit) + IF (explicit) THEN + ! Staging not supported yet, since lowest mode is assumed + ! to be related to centroid + IF (pint_env%transform == transformation_stage) THEN + CPABORT("Centroid constraints are not supported for staging transformation") + END IF + + ! Check thermostats that are not supported: + IF (pint_env%pimd_thermostat == thermostat_gle) THEN + WRITE (UNIT=msg, FMT=*) "GLE Thermostat not supported for "// & + "centroid constraints. Switch to NOSE for numeric integration." + CPABORT(msg) + END IF + ! Warn for NOSE + IF (pint_env%pimd_thermostat == thermostat_nose) THEN + WRITE (UNIT=msg, FMT=*) "PINT WARNING| Nose Thermostat set to "// & + "zero for constrained atoms. Careful interpretation of temperature." + CPWARN(msg) + WRITE (UNIT=msg, FMT=*) "PINT WARNING| Lagrange multipliers are "// & + "are printed every RESPA step and need to be treated carefully." + CPWARN(msg) + END IF + + CALL create_simpar_type(pint_env%simpar) + pint_env%simpar%constraint = explicit + CALL section_vals_val_get(constraint_section, "SHAKE_TOLERANCE", & + r_val=pint_env%simpar%shake_tol) + pint_env%simpar%info_constraint = cp_print_key_unit_nr(pint_env%logger, & + constraint_section, & + "CONSTRAINT_INFO", & + extension=".shakeLog", & + log_filename=.FALSE.) + pint_env%simpar%lagrange_multipliers = cp_print_key_unit_nr(pint_env%logger, & + constraint_section, & + "LAGRANGE_MULTIPLIERS", & + extension=".LagrangeMultLog", & + log_filename=.FALSE.) + pint_env%simpar%dump_lm = & + BTEST(cp_print_key_should_output(pint_env%logger%iter_info, & + constraint_section, & + "LAGRANGE_MULTIPLIERS"), cp_p_file) + + ! Determine constrained atoms: + pint_env%n_atoms_constraints = 0 + DO ig = 1, gci%ncolv%ntot + ! Double counts, if the same atom is involved in different collective variables + pint_env%n_atoms_constraints = pint_env%n_atoms_constraints+SIZE(gci%colv_list(ig)%i_atoms) + END DO + + ALLOCATE (pint_env%atoms_constraints(pint_env%n_atoms_constraints)) + icont = 0 + DO ig = 1, gci%ncolv%ntot + DO iat = 1, SIZE(gci%colv_list(ig)%i_atoms) + icont = icont+1 + pint_env%atoms_constraints(icont) = gci%colv_list(ig)%i_atoms(iat) + END DO + END DO + + ! Set the correction to the temperature due to the frozen degrees of freedom in NOSE: + CALL section_vals_val_get(pint_section, "kT_CORRECTION", & + l_val=ltmp) + IF (ltmp) THEN + pint_env%kTcorr = 1.0_dp+REAL(3*pint_env%n_atoms_constraints, dp)/(REAL(pint_env%ndim, dp)*REAL(pint_env%p, dp)) + ELSE + pint_env%kTcorr = 1.0_dp + END IF + END IF + CALL timestop(handle) RETURN @@ -1125,13 +1216,29 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'pint_init_v', routineP = moduleN//':'//routineN CHARACTER(len=default_string_length) :: msg, stmp, stmp1, stmp2, unit_str - INTEGER :: first_mode, ia, ib, ic, idim, itmp, & - n_rep_val + INTEGER :: first_mode, i, ia, ib, ic, idim, ierr, & + itmp, j, n_rep_val, nparticle, & + nparticle_kind LOGICAL :: done_init, done_quench, done_scale, & done_sped, explicit, ltmp, vels_present - REAL(kind=dp) :: actual_t, ek, rtmp, target_t, unit_conv + REAL(kind=dp) :: actual_t, ek, factor, rtmp, target_t, & + unit_conv + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: vel REAL(kind=dp), DIMENSION(:), POINTER :: r_vals + TYPE(atomic_kind_list_type), POINTER :: atomic_kinds + TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set + TYPE(cell_type), POINTER :: cell TYPE(cp_logger_type), POINTER :: logger + TYPE(cp_subsys_type), POINTER :: subsys + TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles + TYPE(f_env_type), POINTER :: f_env + TYPE(global_constraint_type), POINTER :: gci + TYPE(molecule_kind_list_type), POINTER :: molecule_kinds + TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set + TYPE(molecule_list_type), POINTER :: molecules + TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set + TYPE(particle_list_type), POINTER :: particles + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set TYPE(section_vals_type), POINTER :: input_section CPASSERT(ASSOCIATED(pint_env)) @@ -1140,6 +1247,45 @@ CONTAINS NULLIFY (logger) logger => cp_get_default_logger() + ! Get centroid constraint info, if needed + ! Create a force environment which will be identical to + ! the bead that is being processed by the processor. + IF (pint_env%simpar%constraint) THEN + NULLIFY (subsys, cell) + NULLIFY (atomic_kinds, local_particles, particles) + NULLIFY (local_molecules, molecules, molecule_kinds, gci) + NULLIFY (atomic_kind_set, molecule_kind_set, particle_set, molecule_set) + + CALL f_env_add_defaults(f_env_id=pint_env%replicas%f_env_id, f_env=f_env) + CALL force_env_get(force_env=f_env%force_env, subsys=subsys) + CALL f_env_rm_defaults(f_env, ierr) + CPASSERT(ierr == 0) + + ! Get gci and more from subsys + CALL cp_subsys_get(subsys=subsys, & + cell=cell, & + atomic_kinds=atomic_kinds, & + local_particles=local_particles, & + particles=particles, & + local_molecules=local_molecules, & + molecules=molecules, & + molecule_kinds=molecule_kinds, & + gci=gci) + + nparticle_kind = atomic_kinds%n_els + atomic_kind_set => atomic_kinds%els + molecule_kind_set => molecule_kinds%els + nparticle = particles%n_els + particle_set => particles%els + molecule_set => molecules%els + + ! Allocate work storage + ALLOCATE (vel(3, nparticle)) + vel(:, :) = 0.0_dp + CALL getold(gci, local_molecules, molecule_set, & + molecule_kind_set, particle_set, cell) + END IF + ! read the velocities from the input file if they are given explicitly vels_present = .FALSE. NULLIFY (input_section) @@ -1339,6 +1485,50 @@ CONTAINS END IF END IF + ! Apply constraints to the initial velocities for centroid constraints + IF (pint_env%simpar%constraint) THEN + IF (pint_env%propagator%prop_kind == propagator_rpmd) THEN + ! Multiply with 1/SQRT(n_beads) due to normal mode transformation in RPMD + factor = SQRT(REAL(pint_env%p, dp)) + ELSE + ! lowest NM is centroid + factor = 1 + END IF + + IF (pint_env%logger%para_env%ionode) THEN + DO i = 1, nparticle + DO j = 1, 3 + vel(j, i) = pint_env%uv(1, j+(i-1)*3)/factor + END DO + END DO + + ! Possibly update the target values + CALL shake_update_targets(gci, local_molecules, molecule_set, & + molecule_kind_set, pint_env%dt, & + f_env%force_env%root_section) + CALL rattle_control(gci, local_molecules, molecule_set, & + molecule_kind_set, particle_set, & + vel, pint_env%dt, pint_env%simpar%shake_tol, & + pint_env%simpar%info_constraint, & + pint_env%simpar%lagrange_multipliers, & + .FALSE., & + cell, mp_comm_self, & + local_particles) + END IF + + ! Broadcast updated velocities to other nodes + CALL mp_bcast(vel, & + pint_env%logger%para_env%source, & + pint_env%logger%para_env%group) + + DO i = 1, nparticle + DO j = 1, 3 + pint_env%uv(1, j+(i-1)*3) = vel(j, i)*factor + END DO + END DO + + END IF + RETURN END SUBROUTINE pint_init_v @@ -1523,6 +1713,7 @@ CONTAINS END DO CALL pint_calc_nh_energy(pint_env) END IF + RETURN END SUBROUTINE pint_init_f @@ -1723,6 +1914,7 @@ CONTAINS !> various bug fixes [hforbert] !> 10.2015 Added RPMD propagator and harmonic integrator [Felix Uhl] !> 04.2016 Changed to work with helium_env [cschran] +!> 10.2018 Added centroid constraints [cschran+rperez] !> \author fawzi ! ************************************************************************************************** SUBROUTINE pint_step(pint_env, helium_env) @@ -1732,10 +1924,25 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'pint_step', routineP = moduleN//':'//routineN - INTEGER :: handle, i, ib, idim, inos, iresp - REAL(kind=dp) :: dti, dti2, dti22, e_h, rn, tdti, & - time_start, time_stop, tol + INTEGER :: handle, i, ia, ib, idim, ierr, inos, & + iresp, j, k, nparticle, nparticle_kind + REAL(kind=dp) :: dt_temp, dti, dti2, dti22, e_h, factor, & + rn, tdti, time_start, time_stop, tol + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: pos, vel REAL(kind=dp), DIMENSION(:, :, :), POINTER :: tmp + TYPE(atomic_kind_list_type), POINTER :: atomic_kinds + TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set + TYPE(cell_type), POINTER :: cell + TYPE(cp_subsys_type), POINTER :: subsys + TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles + TYPE(f_env_type), POINTER :: f_env + TYPE(global_constraint_type), POINTER :: gci + TYPE(molecule_kind_list_type), POINTER :: molecule_kinds + TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set + TYPE(molecule_list_type), POINTER :: molecules + TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set + TYPE(particle_list_type), POINTER :: particles + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set CALL timeset(routineN, handle) time_start = m_walltime() @@ -1748,17 +1955,80 @@ CONTAINS tdti = 2.*dti dti22 = dti**2/2._dp + ! Get centroid constraint info, if needed + ! Create a force environment which will be identical to + ! the bead that is being processed by the processor. + IF (pint_env%simpar%constraint) THEN + NULLIFY (subsys, cell) + NULLIFY (atomic_kinds, local_particles, particles) + NULLIFY (local_molecules, molecules, molecule_kinds, gci) + NULLIFY (atomic_kind_set, molecule_kind_set, particle_set, molecule_set) + + CALL f_env_add_defaults(f_env_id=pint_env%replicas%f_env_id, f_env=f_env) + CALL force_env_get(force_env=f_env%force_env, subsys=subsys) + CALL f_env_rm_defaults(f_env, ierr) + CPASSERT(ierr == 0) + + ! Get gci and more from subsys + CALL cp_subsys_get(subsys=subsys, & + cell=cell, & + atomic_kinds=atomic_kinds, & + local_particles=local_particles, & + particles=particles, & + local_molecules=local_molecules, & + molecules=molecules, & + molecule_kinds=molecule_kinds, & + gci=gci) + + nparticle_kind = atomic_kinds%n_els + atomic_kind_set => atomic_kinds%els + molecule_kind_set => molecule_kinds%els + nparticle = particles%n_els + particle_set => particles%els + molecule_set => molecules%els + + ! Allocate work storage + ALLOCATE (pos(3, nparticle)) + ALLOCATE (vel(3, nparticle)) + pos(:, :) = 0.0_dp + vel(:, :) = 0.0_dp + + IF (pint_env%propagator%prop_kind == propagator_rpmd) THEN + ! Multiply with 1/SQRT(n_beads) due to normal mode transformation in RPMD + factor = SQRT(REAL(pint_env%p, dp)) + ELSE + factor = 1 + END IF + + CALL getold(gci, local_molecules, molecule_set, & + molecule_kind_set, particle_set, cell) + END IF + SELECT CASE (pint_env%harm_integrator) CASE (integrate_numeric) + DO iresp = 1, pint_env%nrespa + ! integrate bead positions, first_propagated_mode = { 1, 2 } ! Nose needs an extra step IF (pint_env%pimd_thermostat == thermostat_nose) THEN + + !Set thermostat action of constrained DoF to zero: + IF (pint_env%simpar%constraint) THEN + DO k = 1, pint_env%n_atoms_constraints + ia = pint_env%atoms_constraints(k) + DO j = 3*(ia-1)+1, 3*ia + pint_env%tv(:, 1, j) = 0.0_dp + END DO + END DO + END IF + DO i = pint_env%first_propagated_mode, pint_env%p pint_env%ux(i, :) = pint_env%ux(i, :)- & dti22*pint_env%uv(i, :)*pint_env%tv(1, i, :) END DO pint_env%tx = pint_env%tx+dti*pint_env%tv+dti22*pint_env%tf + END IF !Integrate position in harmonic springs (uf_h) and physical potential !(uf) @@ -1784,12 +2054,64 @@ CONTAINS pint_env%uv_t = pint_env%uv END SELECT - !Integrate hamonic velocities and physical velocities + !Set thermostat action of constrained DoF to zero: + IF (pint_env%simpar%constraint) THEN + DO k = 1, pint_env%n_atoms_constraints + ia = pint_env%atoms_constraints(k) + DO j = 3*(ia-1)+1, 3*ia + pint_env%tv(:, 1, j) = 0.0_dp + pint_env%tv_t(:, 1, j) = 0.0_dp + END DO + END DO + END IF + + !Integrate harmonic velocities and physical velocities pint_env%uv_t = pint_env%uv_t+dti2*(pint_env%uf_h+pint_env%uf) + ! physical forces are only applied in first respa step. pint_env%uf = 0.0_dp ! calc harmonic forces at new pos pint_env%ux = pint_env%ux_t + + ! Apply centroid constraints (SHAKE) + IF (pint_env%simpar%constraint) THEN + IF (pint_env%logger%para_env%ionode) THEN + DO i = 1, nparticle + DO j = 1, 3 + pos(j, i) = pint_env%ux(1, j+(i-1)*3) + vel(j, i) = pint_env%uv_t(1, j+(i-1)*3) + END DO + END DO + + ! Possibly update the target values + CALL shake_update_targets(gci, local_molecules, molecule_set, & + molecule_kind_set, dti, & + f_env%force_env%root_section) + CALL shake_control(gci, local_molecules, molecule_set, & + molecule_kind_set, particle_set, & + pos, vel, dti, pint_env%simpar%shake_tol, & + pint_env%simpar%info_constraint, & + pint_env%simpar%lagrange_multipliers, & + pint_env%simpar%dump_lm, cell, & + mp_comm_self, local_particles) + END IF + ! Positions and velocities of centroid were constrained by SHAKE + CALL mp_bcast(pos, & + pint_env%logger%para_env%source, & + pint_env%logger%para_env%group) + CALL mp_bcast(vel, & + pint_env%logger%para_env%source, & + pint_env%logger%para_env%group) + ! Transform back to normal modes: + DO i = 1, nparticle + DO j = 1, 3 + pint_env%ux(1, j+(i-1)*3) = pos(j, i) + pint_env%uv_t(1, j+(i-1)*3) = vel(j, i) + END DO + END DO + + END IF + CALL pint_calc_uf_h(pint_env=pint_env, e_h=e_h) pint_env%uv_t = pint_env%uv_t+dti2*(pint_env%uf_h+pint_env%uf) @@ -1819,6 +2141,7 @@ CONTAINS !forces pint_env%uf = pint_env%uf*rn pint_env%uv_t = pint_env%uv_t+dti2*pint_env%uf + END IF ! Apply second half of thermostats @@ -1830,11 +2153,21 @@ CONTAINS DO idim = 1, pint_env%ndim DO ib = 1, pint_env%p pint_env%tf(1, ib, idim) = (pint_env%mass_fict(ib, idim)* & - pint_env%uv_new(ib, idim)**2-pint_env%kT)/ & + pint_env%uv_new(ib, idim)**2-pint_env%kT*pint_env%kTcorr)/ & pint_env%Q(ib) END DO END DO + !Set thermostat action of constrained DoF to zero: + IF (pint_env%simpar%constraint) THEN + DO k = 1, pint_env%n_atoms_constraints + ia = pint_env%atoms_constraints(k) + DO j = 3*(ia-1)+1, 3*ia + pint_env%tf(:, 1, j) = 0.0_dp + END DO + END DO + END IF + DO idim = 1, pint_env%ndim DO ib = 1, pint_env%p DO inos = 1, pint_env%nnos-1 @@ -1843,10 +2176,21 @@ CONTAINS (1._dp+dti2*pint_env%tv(inos+1, ib, idim)) pint_env%tf(inos+1, ib, idim) = & (pint_env%tv_new(inos, ib, idim)**2- & - pint_env%kT/pint_env%Q(ib)) + pint_env%kT*pint_env%kTcorr/pint_env%Q(ib)) tol = MAX(tol, ABS(pint_env%tv(inos, ib, idim) & -pint_env%tv_new(inos, ib, idim))) END DO + !Set thermostat action of constrained DoF to zero: + IF (pint_env%simpar%constraint) THEN + DO k = 1, pint_env%n_atoms_constraints + ia = pint_env%atoms_constraints(k) + DO j = 3*(ia-1)+1, 3*ia + pint_env%tv_new(:, 1, j) = 0.0_dp + pint_env%tf(:, 1, j) = 0.0_dp + END DO + END DO + END IF + pint_env%tv_new(pint_env%nnos, ib, idim) = & pint_env%tv_t(pint_env%nnos, ib, idim)+ & dti2*pint_env%tf(pint_env%nnos, ib, idim) @@ -1854,6 +2198,16 @@ CONTAINS -pint_env%tv_new(pint_env%nnos, ib, idim))) tol = MAX(tol, ABS(pint_env%uv(ib, idim) & -pint_env%uv_new(ib, idim))) + !Set thermostat action of constrained DoF to zero: + IF (pint_env%simpar%constraint) THEN + DO k = 1, pint_env%n_atoms_constraints + ia = pint_env%atoms_constraints(k) + DO j = 3*(ia-1)+1, 3*ia + pint_env%tv_new(:, 1, j) = 0.0_dp + END DO + END DO + END IF + END DO END DO @@ -1861,10 +2215,49 @@ CONTAINS pint_env%tv = pint_env%tv_new IF (tol <= pint_env%v_tol) EXIT END DO + + ! Apply centroid constraints (RATTLE) + IF (pint_env%simpar%constraint) THEN + IF (pint_env%logger%para_env%ionode) THEN + DO i = 1, nparticle + DO j = 1, 3 + vel(j, i) = pint_env%uv(1, j+(i-1)*3) + END DO + END DO + + ! Small time step for all small integrations steps + ! Big step for last RESPA + IF (iresp == pint_env%nrespa) THEN + dt_temp = dti + ELSE + dt_temp = dti*rn + END IF + CALL rattle_control(gci, local_molecules, molecule_set, & + molecule_kind_set, particle_set, & + vel, dt_temp, pint_env%simpar%shake_tol, & + pint_env%simpar%info_constraint, & + pint_env%simpar%lagrange_multipliers, & + pint_env%simpar%dump_lm, cell, & + mp_comm_self, local_particles) + END IF + ! Velocities of centroid were constrained by RATTLE + ! Broadcast updated velocities to other nodes + CALL mp_bcast(vel, & + pint_env%logger%para_env%source, & + pint_env%logger%para_env%group) + + DO i = 1, nparticle + DO j = 1, 3 + pint_env%uv(1, j+(i-1)*3) = vel(j, i) + END DO + END DO + END IF + DO inos = 1, pint_env%nnos-1 pint_env%tf(inos, :, :) = pint_env%tf(inos, :, :) & -pint_env%tv(inos, :, :)*pint_env%tv(inos+1, :, :) END DO + CASE (thermostat_gle) CALL pint_gle_step(pint_env) pint_env%uv = pint_env%uv_t @@ -1872,13 +2265,15 @@ CONTAINS pint_env%uv = pint_env%uv_t END SELECT END DO + CASE (integrate_exact) ! The Liouvillian splitting is as follows: ! 1. Thermostat ! 2. 0.5*physical integration - ! 3. Exact harmonic integration + ! 3. Exact harmonic integration ( + apply centroid constraints (SHAKE)) ! 4. 0.5*physical integration ! 5. Thermostat + ! 6. Apply centroid constraints (RATTLE) ! 1. Apply thermostats SELECT CASE (pint_env%pimd_thermostat) @@ -1913,11 +2308,53 @@ CONTAINS ! 3. Exact harmonic integration IF (pint_env%first_propagated_mode == 1) THEN ! The centroid is integrated via standard velocity-verlet + ! Commented out code is only there to show similarities to + ! Numeric integrator pint_env%ux_t(1, :) = pint_env%ux(1, :)+ & - dti*pint_env%uv_t(1, :)+ & - dti22*pint_env%uf_h(1, :) - pint_env%uv_t(1, :) = pint_env%uv_t(1, :)+ & - dti2*pint_env%uf_h(1, :) + dti*pint_env%uv_t(1, :) !+ & + ! dti22*pint_env%uf_h(1, :) + !pint_env%uv_t(1, :) = pint_env%uv_t(1, :)+ & + ! dti2*pint_env%uf_h(1, :) + + ! Apply centroid constraints (SHAKE) + IF (pint_env%simpar%constraint) THEN + IF (pint_env%logger%para_env%ionode) THEN + ! Transform positions and velocities to Cartesian coordinates: + DO i = 1, nparticle + DO j = 1, 3 + pos(j, i) = pint_env%ux_t(1, j+(i-1)*3)/factor + vel(j, i) = pint_env%uv_t(1, j+(i-1)*3)/factor + END DO + END DO + + ! Possibly update the target values + CALL shake_update_targets(gci, local_molecules, molecule_set, & + molecule_kind_set, dti, & + f_env%force_env%root_section) + CALL shake_control(gci, local_molecules, molecule_set, & + molecule_kind_set, particle_set, & + pos, vel, dti, pint_env%simpar%shake_tol, & + pint_env%simpar%info_constraint, & + pint_env%simpar%lagrange_multipliers, & + pint_env%simpar%dump_lm, cell, & + mp_comm_self, local_particles) + END IF + ! Positions and velocities of centroid were constrained by SHAKE + CALL mp_bcast(pos, & + pint_env%logger%para_env%source, & + pint_env%logger%para_env%group) + CALL mp_bcast(vel, & + pint_env%logger%para_env%source, & + pint_env%logger%para_env%group) + ! Transform back to normal modes: + DO i = 1, nparticle + DO j = 1, 3 + pint_env%ux_t(1, j+(i-1)*3) = pos(j, i)*factor + pint_env%uv_t(1, j+(i-1)*3) = vel(j, i)*factor + END DO + END DO + + END IF ELSE ! set velocities to zero for fixed centroids @@ -1936,6 +2373,7 @@ CONTAINS ! 4. 1/2*Physical integration pint_env%uf = 0.0_dp CALL pint_u2x(pint_env) + !TODO apply bead-wise constraints CALL pint_calc_f(pint_env) ! perform helium step and add helium forces IF (PRESENT(helium_env)) THEN @@ -1982,8 +2420,46 @@ CONTAINS pint_env%uv = pint_env%uv_t END SELECT + ! 6. Apply centroid constraints (RATTLE) + IF (pint_env%simpar%constraint) THEN + IF (pint_env%logger%para_env%ionode) THEN + ! Transform positions and velocities to Cartesian coordinates: + DO i = 1, nparticle + DO j = 1, 3 + vel(j, i) = pint_env%uv(1, j+(i-1)*3)/factor + END DO + END DO + + CALL rattle_control(gci, local_molecules, & + molecule_set, molecule_kind_set, & + particle_set, vel, dti, & + pint_env%simpar%shake_tol, & + pint_env%simpar%info_constraint, & + pint_env%simpar%lagrange_multipliers, & + pint_env%simpar%dump_lm, cell, & + mp_comm_self, local_particles) + END IF + ! Velocities of centroid were constrained by RATTLE + ! Broadcast updated velocities to other nodes + CALL mp_bcast(vel, & + pint_env%logger%para_env%source, & + pint_env%logger%para_env%group) + + ! Transform back to normal modes: + ! Multiply with SQRT(n_beads) due to normal mode transformation + DO i = 1, nparticle + DO j = 1, 3 + pint_env%uv(1, j+(i-1)*3) = vel(j, i)*factor + END DO + END DO + + !TODO apply bead-wise constraints + END IF + END SELECT + IF (pint_env%simpar%constraint) DEALLOCATE (pos, vel) + ! calculate the energy components CALL pint_calc_energy(pint_env) CALL pint_calc_total_action(pint_env) diff --git a/src/motion/pint_types.F b/src/motion/pint_types.F index 36c188a6f5..75b0d8c653 100644 --- a/src/motion/pint_types.F +++ b/src/motion/pint_types.F @@ -12,6 +12,7 @@ MODULE pint_types USE parallel_rng_types, ONLY: rng_record_length,& rng_stream_type USE replica_types, ONLY: replica_env_type + USE simpar_types, ONLY: simpar_type #include "../base/base_uses.f90" IMPLICIT NONE @@ -135,6 +136,11 @@ MODULE pint_types TYPE(piglet_therm_type), POINTER :: piglet_therm TYPE(qtb_therm_type), POINTER :: qtb_therm TYPE(pint_propagator_type), POINTER :: propagator + TYPE(simpar_type), POINTER :: simpar + INTEGER :: n_atoms_constraints + INTEGER, DIMENSION(:), POINTER :: atoms_constraints + REAL(KIND=dp) :: kTcorr + END TYPE pint_env_type ! *************************************************************************** @@ -194,13 +200,13 @@ MODULE pint_types !> \author Felix Uhl ! *************************************************************************** TYPE pile_therm_type - INTEGER :: ref_count - REAL(KIND=dp) :: lamb, tau, thermostat_energy - REAL(KIND=dp), DIMENSION(:), POINTER :: c1 - REAL(KIND=dp), DIMENSION(:), POINTER :: c2 - REAL(KIND=dp), DIMENSION(:), POINTER :: g_fric - REAL(KIND=dp), DIMENSION(:, :), POINTER :: massfact - TYPE(rng_stream_type), POINTER :: gaussian_rng_stream + INTEGER :: ref_count + REAL(KIND=dp) :: lamb, tau, thermostat_energy + REAL(KIND=dp), DIMENSION(:), POINTER :: c1 + REAL(KIND=dp), DIMENSION(:), POINTER :: c2 + REAL(KIND=dp), DIMENSION(:), POINTER :: g_fric + REAL(KIND=dp), DIMENSION(:, :), POINTER :: massfact + TYPE(rng_stream_type), POINTER :: gaussian_rng_stream END TYPE pile_therm_type ! *************************************************************************** diff --git a/src/start/input_cp2k_motion.F b/src/start/input_cp2k_motion.F index afa14a8858..121ee94d5e 100644 --- a/src/start/input_cp2k_motion.F +++ b/src/start/input_cp2k_motion.F @@ -1681,6 +1681,12 @@ CONTAINS unit_str="K") CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="kT_CORRECTION", & + description="Corrects for the loss of temperature due to constrained "// & + "degrees of freedom for Nose-Hover chains and numeric integration", & + repeats=.FALSE., default_l_val=.FALSE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) CALL keyword_create(keyword, __LOCATION__, name="T_tol", variants=(/"temp_to"/), & description="threshold for the oscillations of the temperature "// & "excedeed which the temperature is rescaled. 0 means no rescaling.", & diff --git a/tests/Pimd/regtest-1/TEST_FILES b/tests/Pimd/regtest-1/TEST_FILES index 2f9ac0fba3..6e7b8bdc33 100644 --- a/tests/Pimd/regtest-1/TEST_FILES +++ b/tests/Pimd/regtest-1/TEST_FILES @@ -8,6 +8,8 @@ h2o_pint_qs_nose.inp 9 1.0E-14 - h2o_pint_qs_nose_restart.inp 9 2e-09 -17.1295368929579 h2o_pint_exact_harm.inp 9 3e-14 2.8542777100474888E-003 h2o_pint_rpmd.inp 9 2e-14 1.1503487725277912E-002 +h2o_pint_constraints_nose.inp 9 2e-14 9.7527459873828996E-003 +h2o_pint_constraints_exact.inp 9 2e-14 2.8127782657626638E-003 he32_only.inp 40 2e-14 -6.8432379611995943E-005 he32_only_restart.inp 40 2e-14 -6.9355037438372283E-005 water-in-helium.inp 9 7e-14 1.0717546824070511E-003 diff --git a/tests/Pimd/regtest-1/h2o_pint_constraints_exact.inp b/tests/Pimd/regtest-1/h2o_pint_constraints_exact.inp new file mode 100644 index 0000000000..b8ae96fd5f --- /dev/null +++ b/tests/Pimd/regtest-1/h2o_pint_constraints_exact.inp @@ -0,0 +1,67 @@ +&GLOBAL + PROJECT_NAME h2o_pint_constraints_exact + RUN_TYPE PINT + PRINT_LEVEL LOW +&END GLOBAL + +&MOTION + &PINT + P 4 + PROC_PER_REPLICA 1 + NUM_STEPS 4 + DT 0.5 + NRESPA 2 #produces an intended warning. + TEMP 50.0 + TRANSFORMATION NORMAL + HARM_INT EXACT + &END PINT + &CONSTRAINT + &COLLECTIVE + COLVAR 1 + INTERMOLECULAR + &END COLLECTIVE + &LAGRANGE_MULTIPLIERS + COMMON_ITERATION_LEVELS 1 + &END + &END CONSTRAINT +&END MOTION + +&FORCE_EVAL + METHOD FIST + &MM + &FORCEFIELD + PARM_FILE_NAME ../../Fist/sample_pot/water.pot + PARMTYPE CHM + &CHARGE + ATOM OT + CHARGE -0.8476 + &END CHARGE + &CHARGE + ATOM HT + CHARGE 0.4238 + &END CHARGE + &END FORCEFIELD + &POISSON + &EWALD + EWALD_TYPE SPME + ALPHA 0.44 + GMAX 24 + O_SPLINE 6 + &END EWALD + &END POISSON + &END MM + &SUBSYS + &CELL + ABC 24.955 24.955 24.955 + &END CELL + &TOPOLOGY + COORD_FILE_NAME ../../Fist/sample_pdb/water_1.pdb + COORD_FILE_FORMAT PDB + &END TOPOLOGY + &COLVAR + &DISTANCE + ATOMS 1 2 + &END DISTANCE + &END COLVAR + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/Pimd/regtest-1/h2o_pint_constraints_nose.inp b/tests/Pimd/regtest-1/h2o_pint_constraints_nose.inp new file mode 100644 index 0000000000..bb6f2632fe --- /dev/null +++ b/tests/Pimd/regtest-1/h2o_pint_constraints_nose.inp @@ -0,0 +1,69 @@ +&GLOBAL + PROJECT_NAME h2o_pint_constraints_nose + RUN_TYPE PINT + PRINT_LEVEL LOW +&END GLOBAL + +&MOTION + &PINT + P 4 + PROC_PER_REPLICA 1 + NUM_STEPS 4 + DT 0.5 + NRESPA 2 + TEMP 50.0 + TRANSFORMATION NORMAL + &NOSE + NNOS 3 + &END NOSE + &END PINT + &CONSTRAINT + &COLLECTIVE + COLVAR 1 + INTERMOLECULAR + &END COLLECTIVE + &LAGRANGE_MULTIPLIERS + COMMON_ITERATION_LEVELS 1 + &END + &END CONSTRAINT +&END MOTION + +&FORCE_EVAL + METHOD FIST + &MM + &FORCEFIELD + PARM_FILE_NAME ../../Fist/sample_pot/water.pot + PARMTYPE CHM + &CHARGE + ATOM OT + CHARGE -0.8476 + &END CHARGE + &CHARGE + ATOM HT + CHARGE 0.4238 + &END CHARGE + &END FORCEFIELD + &POISSON + &EWALD + EWALD_TYPE SPME + ALPHA 0.44 + GMAX 24 + O_SPLINE 6 + &END EWALD + &END POISSON + &END MM + &SUBSYS + &CELL + ABC 24.955 24.955 24.955 + &END CELL + &TOPOLOGY + COORD_FILE_NAME ../../Fist/sample_pdb/water_1.pdb + COORD_FILE_FORMAT PDB + &END TOPOLOGY + &COLVAR + &DISTANCE + ATOMS 1 2 + &END DISTANCE + &END COLVAR + &END SUBSYS +&END FORCE_EVAL