From ee457522e3acd6ad2eb926a0fac8f467ea757f94 Mon Sep 17 00:00:00 2001 From: Teodoro Laino Date: Tue, 9 Dec 2008 16:18:01 +0000 Subject: [PATCH] VIRTUAL SITE constraint. This allows for TIP4P water. Added regtest (patch created by Marcel Baer and fixed/cleaned by Teo). Reset few regtests for numerics. svn-origin-rev: 8034 --- src/OBJECTDEFS | 1 + src/constraint.F | 46 ++- src/constraint_vsite.F | 293 ++++++++++++++++ src/force_env_methods.F | 4 +- src/input_cp2k_constraints.F | 91 +++++ src/integrator_utils.F | 464 +++++++++++++------------ src/md_vel_utils.F | 62 ++-- src/molecule_kind_types.F | 60 +++- src/molecule_types_new.F | 7 +- src/topology_constraint_util.F | 100 +++++- src/topology_input.F | 94 ++++- src/topology_types.F | 85 +++++ tests/Fist/regtest-10/TEST_FILES_RESET | 4 + tests/Fist/regtest-4/TEST_FILES_RESET | 5 + tests/Fist/regtest-5/TEST_FILES | 2 + tests/Fist/regtest-5/water_tip4p.inp | 224 ++++++++++++ tests/Fist/sample_pdb/water_tip4p.pdb | 10 + tests/Fist/sample_psf/water_tip4p.psf | 41 +++ tests/SE/regtest-2/TEST_FILES_RESET | 2 + 19 files changed, 1307 insertions(+), 288 deletions(-) create mode 100644 src/constraint_vsite.F create mode 100644 tests/Fist/regtest-5/water_tip4p.inp create mode 100644 tests/Fist/sample_pdb/water_tip4p.pdb create mode 100644 tests/Fist/sample_psf/water_tip4p.psf diff --git a/src/OBJECTDEFS b/src/OBJECTDEFS index 90918f84a8..a9df3becae 100644 --- a/src/OBJECTDEFS +++ b/src/OBJECTDEFS @@ -46,6 +46,7 @@ OBJECTS_GENERIC =\ constraint_clv.o\ constraint_fxd.o\ constraint_util.o\ + constraint_vsite.o\ core_ppnl.o\ cp2k_debug.o\ cp2k_info.o\ diff --git a/src/constraint.F b/src/constraint.F index f8b5bb5906..a536704855 100644 --- a/src/constraint.F +++ b/src/constraint.F @@ -37,6 +37,8 @@ MODULE constraint get_roll_matrix,& restore_temporary_set,& update_temporary_set + USE constraint_vsite, ONLY: shake_vsite_ext,& + shake_vsite_int USE cp_para_types, ONLY: cp_para_env_type USE distribution_1d_types, ONLY: distribution_1d_type USE f77_blas @@ -100,7 +102,7 @@ CONTAINS routineP = moduleN//':'//routineN INTEGER :: handle, i, ikind, imol, ishake_ext, ishake_int, k, n3x3con, & - n4x6con, nconstraint, nkind, nmol_per_kind + n4x6con, nconstraint, nkind, nmol_per_kind, nvsitecon LOGICAL :: do_ext_constraint REAL(KIND=dp) :: int_max_sigma, mass, max_sigma REAL(KIND=dp), DIMENSION(SIZE(pos, 2)) :: imass @@ -114,7 +116,7 @@ CONTAINS DO k = 1, SIZE(pos,2) atomic_kind=>particle_set (k) % atomic_kind CALL get_atomic_kind(atomic_kind=atomic_kind,mass=mass) - imass(k) = 1.0_dp/mass + imass(k) = 1.0_dp/MAX(mass,EPSILON(0.0_dp)) END DO do_ext_constraint = (gci%ntot/=0) ishake_ext = 0 @@ -131,7 +133,7 @@ CONTAINS molecule_kind => molecule % molecule_kind CALL get_molecule_kind ( molecule_kind, ncolv = ncolv, & ng3x3 = n3x3con, ng4x6=n4x6con,& - nconstraint=nconstraint ) + nconstraint=nconstraint,nvsite=nvsitecon) IF (nconstraint==0) CYCLE ishake_int = 0 int_max_sigma = -1.0E+10_dp @@ -153,6 +155,10 @@ CONTAINS END DO Shake_Intra_Loop max_sigma = MAX(max_sigma, int_max_sigma) CALL shake_int_info(log_unit,i,ishake_int,max_sigma,error) + ! Virtual Site + IF ( nvsitecon /= 0 ) & + CALL shake_vsite_int( molecule, particle_set, pos, vel, dt, ishake_int,& + int_max_sigma, error) END DO END DO MOL ! Intermolecular constraints @@ -170,12 +176,17 @@ CONTAINS IF ( gci%ncolv%ntot /= 0 ) & CALL shake_colv_ext(gci, particle_set, pos, vel, dt, ishake_ext,& cell, imass, max_sigma, error) + ! Virtual Site + IF ( gci%nvsite /= 0 ) & + CALL shake_vsite_ext( gci, particle_set, pos, vel, dt, ishake_ext,& + max_sigma, error) CALL restore_temporary_set( particle_set, local_particles, pos=pos, vel=vel ) END IF CALL shake_ext_info(log_unit,ishake_ext,max_sigma,error) END DO Shake_Inter_Loop CALL dump_lagrange_mult(dump_lm, lagrange_mult,local_molecules,molecule_set,gci,& molecule_kind_set,group,"S",error) + CALL timestop(handle) END SUBROUTINE shake_control @@ -219,7 +230,7 @@ CONTAINS DO k = 1, SIZE(vel,2) atomic_kind=>particle_set (k) % atomic_kind CALL get_atomic_kind(atomic_kind=atomic_kind,mass=mass) - imass(k) = 1.0_dp/mass + imass(k) = 1.0_dp/MAX(mass,EPSILON(0.0_dp)) END DO do_ext_constraint = (gci%ntot/=0) irattle_ext = 0 @@ -311,7 +322,7 @@ CONTAINS INTEGER :: handle, i, ikind, imol, ishake_ext, ishake_int, k, & lagrange_mult, log_unit, n3x3con, n4x6con, nconstraint, nkind, & - nmol_per_kind + nmol_per_kind, nvsitecon LOGICAL :: do_ext_constraint, dump_lm REAL(KIND=dp) :: int_max_sigma, mass, & max_sigma, shake_tol @@ -337,7 +348,7 @@ CONTAINS DO k = 1, SIZE(pos,2) atomic_kind=>particle_set (k) % atomic_kind CALL get_atomic_kind(atomic_kind=atomic_kind,mass=mass) - imass(k) = 1.0_dp/mass + imass(k) = 1.0_dp/MAX(mass,EPSILON(0.0_dp)) END DO do_ext_constraint = (gci%ntot/=0) ishake_ext = 0 @@ -354,7 +365,7 @@ CONTAINS molecule_kind => molecule % molecule_kind CALL get_molecule_kind ( molecule_kind, ncolv = ncolv, & ng3x3 = n3x3con, ng4x6=n4x6con,& - nconstraint=nconstraint ) + nconstraint=nconstraint,nvsite=nvsitecon) IF (nconstraint==0) CYCLE ishake_int = 0 int_max_sigma = -1.0E+10_dp @@ -376,6 +387,12 @@ CONTAINS END DO Shake_Roll_Intra_Loop max_sigma = MAX(max_sigma, int_max_sigma) CALL shake_int_info(log_unit,i,ishake_int,max_sigma,error) + ! Virtual Site + IF ( nvsitecon /= 0 ) THEN + CALL cp_unimplemented_error(fromWhere=routineP, & + message="Virtual Site Constraint/Restraint not implemented for SHAKE_ROLL!", & + error=error, error_level=cp_failure_level) + END IF END DO END DO MOL ! Intermolecular constraints @@ -393,6 +410,11 @@ CONTAINS IF ( gci%ncolv%ntot /= 0 ) & CALL shake_roll_colv_ext( gci, particle_set, pos, vel, r_shake,& v_shake, dt, ishake_ext, cell, imass, max_sigma, error ) + ! Virtual Site + IF ( gci%nvsite /= 0 ) & + CALL cp_unimplemented_error(fromWhere=routineP, & + message="Virtual Site Constraint/Restraint not implemented for SHAKE_ROLL!", & + error=error, error_level=cp_failure_level) CALL restore_temporary_set ( particle_set, local_particles, pos = pos, vel = vel ) END IF CALL shake_ext_info(log_unit,ishake_ext,max_sigma,error) @@ -410,7 +432,7 @@ CONTAINS ! ***************************************************************************** SUBROUTINE rattle_roll_control( gci, local_molecules, molecule_set, & molecule_kind_set, particle_set, vel, dt, simpar, virial, vector,& - veps, mass, roll_tol, iroll, para_env, u, cell, local_particles, error ) + veps, roll_tol, iroll, para_env, u, cell, local_particles, error ) TYPE(global_constraint_type), POINTER :: gci TYPE(distribution_1d_type), POINTER :: local_molecules @@ -424,8 +446,6 @@ CONTAINS REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: vector REAL(KIND=dp), DIMENSION(:, :), & INTENT(INOUT) :: veps - REAL(KIND=dp), DIMENSION(:, :), & - INTENT(IN) :: mass REAL(KIND=dp), INTENT(OUT) :: roll_tol INTEGER, INTENT(INOUT) :: iroll TYPE(cp_para_env_type), INTENT(IN) :: para_env @@ -442,7 +462,7 @@ CONTAINS lagrange_mult, log_unit, n3x3con, n4x6con, nconstraint, nkind, & nmol_per_kind LOGICAL :: do_ext_constraint, dump_lm - REAL(KIND=dp) :: amass, int_max_sigma, & + REAL(KIND=dp) :: int_max_sigma, mass, & max_sigma, rattle_tol REAL(KIND=dp), DIMENSION(3, 3) :: r_rattle REAL(KIND=dp), DIMENSION(SIZE(vel, 2)) :: imass @@ -466,8 +486,8 @@ CONTAINS END IF DO k = 1, SIZE(vel,2) atomic_kind=>particle_set (k) % atomic_kind - CALL get_atomic_kind(atomic_kind=atomic_kind,mass=amass) - imass(k) = 1.0_dp/amass + CALL get_atomic_kind(atomic_kind=atomic_kind,mass=mass) + imass(k) = 1.0_dp/MAX(mass,EPSILON(0.0_dp)) END DO do_ext_constraint = (gci%ntot/=0) irattle_ext = 0 diff --git a/src/constraint_vsite.F b/src/constraint_vsite.F new file mode 100644 index 0000000000..2b48c6a9fb --- /dev/null +++ b/src/constraint_vsite.F @@ -0,0 +1,293 @@ +!-----------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright (C) 2000 - 2008 CP2K developers group ! +!-----------------------------------------------------------------------------! + +! ***************************************************************************** +!> \brief Routines to handle the virtual site constraint/restraint +!> \par History +!> Teodoro Laino [tlaino] 12.2008 - Preparing for VIRTUAL SITE constraints +!> (patch by Marcel Baer) +! ***************************************************************************** +MODULE constraint_vsite + USE cp_subsys_types, ONLY: cp_subsys_get,& + cp_subsys_type + USE distribution_1d_types, ONLY: distribution_1d_type + USE f77_blas + USE force_env_types, ONLY: force_env_get,& + force_env_type + USE kinds, ONLY: dp + USE mol_kind_new_list_types, ONLY: mol_kind_new_list_type + USE mol_new_list_types, ONLY: mol_new_list_type + USE molecule_kind_types, ONLY: fixd_constraint_type,& + get_molecule_kind,& + molecule_kind_type,& + vsite_constraint_type + USE molecule_types_new, ONLY: get_molecule,& + global_constraint_type,& + molecule_type + USE particle_list_types, ONLY: particle_list_type + USE particle_types, ONLY: particle_type +#include "cp_common_uses.h" + + IMPLICIT NONE + + PRIVATE + PUBLIC :: shake_vsite_int,& + force_vsite_int,& + force_vsite_ext,& + shake_vsite_ext,& + vsite_force_control + + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'constraint_vsite' + +CONTAINS + +! ***************************************************************************** +!> \brief control force distribution for virtual sites +!> \par History +!> - none +!> \author Marcel Baer +!> \date 12.2008 +! ***************************************************************************** + SUBROUTINE vsite_force_control( force_env, error,w) + TYPE(force_env_type), POINTER :: force_env + TYPE(cp_error_type), INTENT(inout) :: error + REAL(KIND=dp), DIMENSION(:, :), OPTIONAL :: w + + INTEGER :: i, ikind, imol, nconstraint, & + nkind, nmol_per_kind, & + nvsitecon + LOGICAL :: do_ext_constraint + TYPE(cp_subsys_type), POINTER :: subsys + TYPE(distribution_1d_type), POINTER :: local_molecules, & + local_particles + TYPE(global_constraint_type), POINTER :: gci + TYPE(mol_kind_new_list_type), POINTER :: molecule_kinds + TYPE(mol_new_list_type), POINTER :: molecules + TYPE(molecule_kind_type), DIMENSION(:), & + POINTER :: molecule_kind_set + TYPE(molecule_kind_type), POINTER :: molecule_kind + TYPE(molecule_type), DIMENSION(:), & + POINTER :: molecule_set + TYPE(molecule_type), POINTER :: molecule + TYPE(particle_list_type), POINTER :: particles + TYPE(particle_type), DIMENSION(:), & + POINTER :: particle_set + + NULLIFY(gci, subsys, local_molecules, local_particles,& + molecule_kinds) + + CALL force_env_get(force_env=force_env, subsys=subsys,error=error) + + CALL cp_subsys_get(subsys=subsys,local_particles=local_particles,& + particles=particles,local_molecules_new=local_molecules,& + molecule_kinds_new=molecule_kinds,gci=gci,molecules_new=molecules,& + error=error) + + molecule_kind_set => molecule_kinds%els + molecule_set => molecules%els + particle_set => particles%els + nkind = SIZE ( molecule_kind_set ) + ! Intermolecular Virtual Site Constraints + do_ext_constraint = .FALSE. + IF (ASSOCIATED(gci)) THEN + do_ext_constraint = (gci%ntot/=0) + END IF + ! Intramolecular Virtual Site Constraints + MOL: DO ikind = 1, nkind + nmol_per_kind = local_molecules % n_el ( ikind ) + DO imol = 1, nmol_per_kind + i = local_molecules % list ( ikind ) % array ( imol ) + molecule => molecule_set ( i ) + molecule_kind => molecule % molecule_kind + CALL get_molecule_kind ( molecule_kind,nconstraint=nconstraint,nvsite=nvsitecon) + IF (nconstraint==0) CYCLE + IF (nvsitecon/=0) & + CALL force_vsite_int( molecule, particle_set) + END DO + END DO MOL + ! Intermolecular Virtual Site Constraints + IF (do_ext_constraint) THEN + IF ( gci%nvsite /= 0 ) & + CALL force_vsite_ext( gci, particle_set) + END IF + + END SUBROUTINE vsite_force_control + +! ***************************************************************************** +!> \brief Intramolecular virtual site +!> \par History +!> 12.2008 Marcel Baer +! ***************************************************************************** + SUBROUTINE shake_vsite_int( molecule, particle_set, pos, vel, dt, ishake,& + max_sigma, error) + TYPE(molecule_type), POINTER :: molecule + TYPE(particle_type), POINTER :: particle_set( : ) + REAL(KIND=dp), INTENT(INOUT) :: pos( :, : ), vel( :, : ) + REAL(kind=dp), INTENT(in) :: dt + INTEGER, INTENT(IN) :: ishake + REAL(KIND=dp), INTENT(INOUT) :: max_sigma + TYPE(cp_error_type), INTENT(inout) :: error + + INTEGER :: first_atom, nvsite + TYPE(fixd_constraint_type), & + DIMENSION(:), POINTER :: fixd_list + TYPE(molecule_kind_type), POINTER :: molecule_kind + TYPE(vsite_constraint_type), POINTER :: vsite_list( : ) + + NULLIFY(fixd_list) + molecule_kind => molecule % molecule_kind + CALL get_molecule_kind ( molecule_kind, nvsite = nvsite,& + fixd_list=fixd_list, vsite_list = vsite_list ) + CALL get_molecule ( molecule, first_atom = first_atom) + ! Real Shake + CALL shake_vsite_low( fixd_list, vsite_list, nvsite, first_atom, & + particle_set, pos, vel, dt, ishake, max_sigma, error) + + END SUBROUTINE shake_vsite_int + +! ***************************************************************************** +!> \brief Intramolecular virtual site +!> \par History +!> 12.2008 Marcel Baer +! ***************************************************************************** + SUBROUTINE shake_vsite_ext( gci, particle_set, pos, vel, dt, ishake,& + max_sigma, error) + + TYPE(global_constraint_type), POINTER :: gci + TYPE(particle_type), POINTER :: particle_set( : ) + REAL(KIND=dp), INTENT(INOUT) :: pos( :, : ), vel( :, : ) + REAL(kind=dp), INTENT(in) :: dt + INTEGER, INTENT(IN) :: ishake + REAL(KIND=dp), INTENT(INOUT) :: max_sigma + TYPE(cp_error_type), INTENT(inout) :: error + + INTEGER :: first_atom, nvsite + TYPE(fixd_constraint_type), & + DIMENSION(:), POINTER :: fixd_list + TYPE(vsite_constraint_type), POINTER :: vsite_list( : ) + + first_atom = 1 + nvsite = gci%nvsite + fixd_list => gci%fixd_list + vsite_list => gci%vsite_list + ! Real Shake + CALL shake_vsite_low( fixd_list, vsite_list, nvsite, first_atom, & + particle_set, pos, vel, dt, ishake, max_sigma, error) + + END SUBROUTINE shake_vsite_ext + +! ***************************************************************************** +!> \par History +!> 12.2008 Marcel Bear +! ***************************************************************************** + SUBROUTINE shake_vsite_low(fixd_list, vsite_list, nvsite, first_atom, & + particle_set, pos, vel, dt, ishake, max_sigma, error) + TYPE(fixd_constraint_type), & + DIMENSION(:), POINTER :: fixd_list + TYPE(vsite_constraint_type), POINTER :: vsite_list( : ) + INTEGER, INTENT(IN) :: nvsite, first_atom + TYPE(particle_type), POINTER :: particle_set( : ) + REAL(KIND=dp), INTENT(INOUT) :: pos( :, : ), vel( :, : ) + REAL(kind=dp), INTENT(in) :: dt + INTEGER, INTENT(IN) :: ishake + REAL(KIND=dp), INTENT(INOUT) :: max_sigma + TYPE(cp_error_type), INTENT(inout) :: error + + INTEGER :: iconst, index_a, index_b, & + index_c, index_d + REAL(KIND=dp), DIMENSION(3) :: r1, r2 + + DO iconst = 1, nvsite + IF (vsite_list(iconst)%restraint%active) CYCLE + index_a = vsite_list(iconst)%a + first_atom -1 + index_b = vsite_list(iconst)%b + first_atom -1 + index_c = vsite_list(iconst)%c + first_atom -1 + index_d = vsite_list(iconst)%d + first_atom -1 + + r1(:) = pos(:,index_b) - pos(:,index_c ) + r2(:) = pos(:,index_d) - pos(:,index_c ) + pos(:,index_a)=pos(:,index_c)+vsite_list(iconst)%wbc*r1(:)+& + vsite_list(iconst)%wdc*r2(:) + END DO + END SUBROUTINE shake_vsite_low + +! ***************************************************************************** +!> \brief Intramolecular virtual site +!> \par History +!> 12.2008 Marcel Bear +! ***************************************************************************** + SUBROUTINE force_vsite_int( molecule, particle_set) + TYPE(molecule_type), POINTER :: molecule + TYPE(particle_type), POINTER :: particle_set( : ) + + INTEGER :: first_atom, iconst, index_a, & + index_b, index_c, index_d, & + nvsite + REAL(KIND=dp) :: wb, wc, wd + TYPE(molecule_kind_type), POINTER :: molecule_kind + TYPE(vsite_constraint_type), POINTER :: vsite_list( : ) + + molecule_kind => molecule % molecule_kind + CALL get_molecule_kind ( molecule_kind, nvsite = nvsite, vsite_list = vsite_list ) + CALL get_molecule ( molecule, first_atom = first_atom) + + DO iconst = 1, nvsite + IF (vsite_list(iconst)%restraint%active) CYCLE + index_a = vsite_list(iconst)%a + first_atom -1 + index_b = vsite_list(iconst)%b + first_atom -1 + index_c = vsite_list(iconst)%c + first_atom -1 + index_d = vsite_list(iconst)%d + first_atom -1 + + wb=vsite_list(iconst)%wbc + wd=vsite_list(iconst)%wdc + wc=1.0_dp-vsite_list(iconst)%wbc-vsite_list(iconst)%wdc + + particle_set(index_b)%f(:)=particle_set(index_b)%f(:) + wb*particle_set(index_a)%f(:) + particle_set(index_c)%f(:)=particle_set(index_c)%f(:) + wc*particle_set(index_a)%f(:) + particle_set(index_d)%f(:)=particle_set(index_d)%f(:) + wd*particle_set(index_a)%f(:) + particle_set(index_a)%f(:)=0.0_dp + END DO + + END SUBROUTINE force_vsite_int + +! ***************************************************************************** +!> \brief Intramolecular virtual site +!> \par History +!> 12.2008 Marcel Bear +! ***************************************************************************** + SUBROUTINE force_vsite_ext( gci, particle_set) + TYPE(global_constraint_type), POINTER :: gci + TYPE(particle_type), POINTER :: particle_set( : ) + + INTEGER :: first_atom, iconst, index_a, & + index_b, index_c, index_d, & + nvsite + REAL(KIND=dp) :: wb, wc, wd + TYPE(vsite_constraint_type), POINTER :: vsite_list( : ) + + first_atom = 1 + nvsite = gci%nvsite + vsite_list => gci%vsite_list + ! Real Shake + + DO iconst = 1, nvsite + IF (vsite_list(iconst)%restraint%active) CYCLE + index_a = vsite_list(iconst)%a + first_atom -1 + index_b = vsite_list(iconst)%b + first_atom -1 + index_c = vsite_list(iconst)%c + first_atom -1 + index_d = vsite_list(iconst)%d + first_atom -1 + + wb=vsite_list(iconst)%wbc + wd=vsite_list(iconst)%wdc + wc=1.0_dp-vsite_list(iconst)%wbc-vsite_list(iconst)%wdc + + particle_set(index_b)%f(:)=particle_set(index_b)%f(:) + wb*particle_set(index_a)%f(:) + particle_set(index_c)%f(:)=particle_set(index_c)%f(:) + wc*particle_set(index_a)%f(:) + particle_set(index_d)%f(:)=particle_set(index_d)%f(:) + wd*particle_set(index_a)%f(:) + particle_set(index_a)%f(:)=0.0_dp + END DO + END SUBROUTINE force_vsite_ext + +END MODULE constraint_vsite diff --git a/src/force_env_methods.F b/src/force_env_methods.F index 15a27723e8..cc4b41f262 100644 --- a/src/force_env_methods.F +++ b/src/force_env_methods.F @@ -20,6 +20,7 @@ MODULE force_env_methods shake_control USE constraint_fxd, ONLY: fix_atom_control USE constraint_util, ONLY: getold + USE constraint_vsite, ONLY: vsite_force_control USE cp_output_handling, ONLY: cp_p_file,& cp_print_key_finished_output,& cp_print_key_should_output,& @@ -299,6 +300,8 @@ CONTAINS CALL fix_atom_control(force_env, error=error) ! All Restraints CALL restraint_control(force_env, error=error) + ! Virtual Sites + CALL vsite_force_control(force_env,error) ! External Potential CALL add_external_potential(force_env, error=error) ! Rescale forces if requested @@ -1488,7 +1491,6 @@ CPSourceFileRef,& TYPE(particle_list_type), POINTER :: particles failure=.FALSE. - CALL timeset(routineN,handle) CPPrecondition(ASSOCIATED(force_env),cp_failure_level,routineP,error,failure) CPPrecondition(force_env%ref_count>0,cp_failure_level,routineP,error,failure) diff --git a/src/input_cp2k_constraints.F b/src/input_cp2k_constraints.F index 41c4c87156..ca3216d9a8 100644 --- a/src/input_cp2k_constraints.F +++ b/src/input_cp2k_constraints.F @@ -6,6 +6,8 @@ ! ***************************************************************************** !> \par History !> 10.2005 split input_cp2k into smaller modules [fawzi] +!> Teodoro Laino [tlaino] 12.2008 - Preparing for VIRTUAL SITE constraints +!> (patch by Marcel Baer) !> \author teo & fawzi ! ***************************************************************************** MODULE input_cp2k_constraints @@ -106,6 +108,11 @@ CONTAINS CALL restraint_info_section(subsection, error=error) CALL section_add_subsection(section, subsection, error=error) CALL section_release(subsection,error=error) + + CALL create_vsite_section(subsection, error=error) + CALL restraint_info_section(subsection, error=error) + CALL section_add_subsection(section, subsection, error=error) + CALL section_release(subsection,error=error) CALL create_collective_section(subsection, error=error) CALL restraint_info_section(subsection, error=error) @@ -708,4 +715,88 @@ CONTAINS END SUBROUTINE create_g4x6_section +! ***************************************************************************** +!> \brief Create the constraint section specialized on vsite constraints +!> \param section the section to create +!> \param error variable to control error logging, stopping,... +!> see module cp_error_handling +!> \author marcel baer +! ***************************************************************************** + SUBROUTINE create_vsite_section(section,error) + TYPE(section_type), POINTER :: section + TYPE(cp_error_type), INTENT(inout) :: error + + CHARACTER(len=*), PARAMETER :: routineN = 'create_vsite_section', & + routineP = moduleN//':'//routineN + + LOGICAL :: failure + TYPE(keyword_type), POINTER :: keyword + + failure=.FALSE. + + CPPrecondition(.NOT.ASSOCIATED(section),cp_failure_level,routineP,error,failure) + IF (.NOT. failure) THEN + CALL section_create(section,name="virtual_site",& + description="This section is used to set a virtual interaction-site constraint.",& + n_keywords=3, n_subsections=0, repeats=.TRUE., required=.TRUE.,& + error=error) + + NULLIFY(keyword) + + ! Integer + CALL keyword_create(keyword, name="MOLECULE",& + variants=(/"MOL"/),& + description="Specifies the molecule number on which constraint will be applied."//& + " MOLECULE and MOLNAME keyword exclude themself mutually.",& + usage="MOL {integer}", required=.TRUE., n_var=1, type_of_var=integer_t, error=error) + CALL section_add_keyword(section,keyword,error=error) + CALL keyword_release(keyword,error=error) + + CALL keyword_create(keyword, name="MOLNAME",& + variants=(/"SEGNAME"/),& + description="Specifies the name of the molecule on which the constraint will be applied.",& + usage="MOLNAME {character}", required=.TRUE., n_var=1, type_of_var=char_t, error=error) + CALL section_add_keyword(section,keyword,error=error) + CALL keyword_release(keyword,error=error) + + CALL keyword_create(keyword, name="INTERMOLECULAR",& + description="Specify if the constraint/restraint is intermolecular.",& + usage="INTERMOLECULAR ",& + default_l_val=.FALSE.,lone_keyword_l_val=.TRUE.,error=error) + CALL section_add_keyword(section,keyword,error=error) + CALL keyword_release(keyword,error=error) + + CALL keyword_create(keyword, name="ATOMS",& + description="Atoms' index on which apply the constraint (v i j k), first is virtual site",& + usage="ATOMS 1 2 3 4",& + n_var=4,type_of_var=integer_t,error=error) + CALL section_add_keyword(section,keyword,error=error) + CALL keyword_release(keyword,error=error) + + ! Real + CALL keyword_create(keyword, name="PARAMETERS",& + description="The constrained paramters' values to construct virtual site."//& + "r_v=a*r_ij+b*r_kj",& + usage="PARAMETERS {real} {real}",& + type_of_var=real_t, n_var=2,unit_str="internal_cp2k",error=error) + CALL section_add_keyword(section,keyword,error=error) + CALL keyword_release(keyword,error=error) + + ! Logical + CALL keyword_create(keyword, name="EXCLUDE_QM",& + description="Does not apply the constraint to the QM region within a QM/MM calculation",& + usage="EXCLUDE_QM ",& + default_l_val=.FALSE.,lone_keyword_l_val=.TRUE.,error=error) + CALL section_add_keyword(section,keyword,error=error) + CALL keyword_release(keyword,error=error) + + CALL keyword_create(keyword, name="EXCLUDE_MM",& + description="Does not apply the constraint to the MM region within a QM/MM calculation",& + usage="EXCLUDE_MM ",& + default_l_val=.FALSE.,lone_keyword_l_val=.TRUE.,error=error) + CALL section_add_keyword(section,keyword,error=error) + CALL keyword_release(keyword,error=error) + END IF + + END SUBROUTINE create_vsite_section END MODULE input_cp2k_constraints diff --git a/src/integrator_utils.F b/src/integrator_utils.F index 6845e4dd94..e11ba842ab 100644 --- a/src/integrator_utils.F +++ b/src/integrator_utils.F @@ -7,6 +7,8 @@ !> \brief Provides integrator utility routines for the integrators !> \par History !> Teodoro Laino [tlaino] 02.2008 - Splitting from integrator and creation +!> Teodoro Laino [tlaino] 12.2008 - Preparing for VIRTUAL SITE constraints +!> (patch by Marcel Baer) ! ***************************************************************************** MODULE integrator_utils @@ -375,21 +377,23 @@ CONTAINS DO iparticle_kind=1,nparticle_kind atomic_kind => atomic_kind_set(iparticle_kind) CALL get_atomic_kind(atomic_kind=atomic_kind,mass=mass) - nparticle_local = local_particles%n_el(iparticle_kind) - IF(my_tmpv) THEN - DO iparticle_local=1,nparticle_local - iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) - K=K+ 0.5_dp*mass *DOT_PRODUCT(tmp%vel(:,iparticle),tmp%vel(:,iparticle)) - a=a+ DOT_PRODUCT(tmp%vel(:,iparticle),particle_set(iparticle)%f(:)) - b=b+(1.0_dp/mass)*DOT_PRODUCT(particle_set(iparticle)%f(:),particle_set(iparticle)%f(:)) - END DO - ELSE - DO iparticle_local=1,nparticle_local - iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) - K=K+ 0.5_dp*mass *DOT_PRODUCT(particle_set(iparticle)%v(:),particle_set(iparticle)%v(:)) - a=a+ DOT_PRODUCT(particle_set(iparticle)%v(:),particle_set(iparticle)%f(:)) - b=b+(1.0_dp/mass)*DOT_PRODUCT(particle_set(iparticle)%f(:),particle_set(iparticle)%f(:)) - END DO + IF (mass/=0.0_dp) THEN + nparticle_local = local_particles%n_el(iparticle_kind) + IF(my_tmpv) THEN + DO iparticle_local=1,nparticle_local + iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) + K=K+ 0.5_dp*mass *DOT_PRODUCT(tmp%vel(:,iparticle),tmp%vel(:,iparticle)) + a=a+ DOT_PRODUCT(tmp%vel(:,iparticle),particle_set(iparticle)%f(:)) + b=b+(1.0_dp/mass)*DOT_PRODUCT(particle_set(iparticle)%f(:),particle_set(iparticle)%f(:)) + END DO + ELSE + DO iparticle_local=1,nparticle_local + iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) + K=K+ 0.5_dp*mass *DOT_PRODUCT(particle_set(iparticle)%v(:),particle_set(iparticle)%v(:)) + a=a+ DOT_PRODUCT(particle_set(iparticle)%v(:),particle_set(iparticle)%f(:)) + b=b+(1.0_dp/mass)*DOT_PRODUCT(particle_set(iparticle)%f(:),particle_set(iparticle)%f(:)) + END DO + END IF END IF END DO CALL mp_sum(K,para_env % group) @@ -746,8 +750,7 @@ CONTAINS ! resetting CALL set ( old, atomic_kind_set, particle_set, vel, local_particles, cell , npt, 'B' ) CALL rattle_roll_control( gci, local_molecules, molecule_set, molecule_kind_set, & - particle_set, vel, dt, simpar, virial, vector_v, & - npt_loc % v, npt_loc % mass, roll_tol, iroll, & + particle_set, vel, dt, simpar, virial, vector_v, npt_loc%v, roll_tol, iroll,& para_env, u, cell, local_particles, error=error ) END SUBROUTINE rattle_roll_setup @@ -1011,140 +1014,154 @@ CONTAINS DO iparticle_kind=1,nparticle_kind atomic_kind => atomic_kind_set(iparticle_kind) CALL get_atomic_kind(atomic_kind=atomic_kind,mass=mass,shell_active=is_shell) - dm = 0.5_dp*dt/mass - IF(is_shell .AND. shell_adiabatic) THEN - CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell) - dms = 0.5_dp*dt/shell%mass_shell - dmc = 0.5_dp*dt/shell%mass_core - fac_masss = shell%mass_shell/mass - fac_massc = shell%mass_core/mass - nparticle_local = local_particles%n_el(iparticle_kind) - IF(PRESENT(u)) THEN - DO iparticle_local=1,nparticle_local - iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) - shell_index = particle_set(iparticle)%shell_index - ! Transform positions and velocities and forces of the shell - CALL transform_first(shell_particle_set,tmp%shell_pos,tmp%shell_vel,& - shell_index,u,dms,dt,tmp%poly_v,tmp%poly_r,tmp%scale_v,tmp%scale_r) + IF (mass/=0.0_dp) THEN + dm = 0.5_dp*dt/mass + IF(is_shell .AND. shell_adiabatic) THEN + CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell) + dms = 0.5_dp*dt/shell%mass_shell + dmc = 0.5_dp*dt/shell%mass_core + fac_masss = shell%mass_shell/mass + fac_massc = shell%mass_core/mass + nparticle_local = local_particles%n_el(iparticle_kind) + IF(PRESENT(u)) THEN + DO iparticle_local=1,nparticle_local + iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) + shell_index = particle_set(iparticle)%shell_index + ! Transform positions and velocities and forces of the shell + CALL transform_first(shell_particle_set,tmp%shell_pos,tmp%shell_vel,& + shell_index,u,dms,dt,tmp%poly_v,tmp%poly_r,tmp%scale_v,tmp%scale_r) - ! Transform positions and velocities and forces of the core - CALL transform_first(core_particle_set,tmp%core_pos,tmp%core_vel,& - shell_index,u,dmc,dt,tmp%poly_v,tmp%poly_r,tmp%scale_v,tmp%scale_r) + ! Transform positions and velocities and forces of the core + CALL transform_first(core_particle_set,tmp%core_pos,tmp%core_vel,& + shell_index,u,dmc,dt,tmp%poly_v,tmp%poly_r,tmp%scale_v,tmp%scale_r) - ! Derive velocities and positions of the COM - tmp%vel (:,iparticle) = fac_masss*tmp%shell_vel (:,shell_index)+& - fac_massc*tmp%core_vel (:,shell_index) - tmp%pos (:,iparticle) = fac_masss*tmp%shell_pos (:,shell_index)+& - fac_massc*tmp%core_pos (:,shell_index) + ! Derive velocities and positions of the COM + tmp%vel (:,iparticle) = fac_masss*tmp%shell_vel (:,shell_index)+& + fac_massc*tmp%core_vel (:,shell_index) + tmp%pos (:,iparticle) = fac_masss*tmp%shell_pos (:,shell_index)+& + fac_massc*tmp%core_pos (:,shell_index) - tmp%max_vel = MAX(tmp%max_vel,ABS(tmp%vel(1,iparticle)),ABS(tmp%vel (2,iparticle)),ABS(tmp%vel (3,iparticle))) - tmp%max_vel_sc = MAX(tmp%max_vel_sc,ABS(tmp%shell_vel(1,shell_index)-tmp%core_vel(1,shell_index)),& - ABS(tmp%shell_vel(2,shell_index)-tmp%core_vel(2,shell_index)),& - ABS(tmp%shell_vel(3,shell_index)-tmp%core_vel(3,shell_index))) - tmp%max_dr = MAX(tmp%max_dr,ABS(particle_set(iparticle)%r(1)-tmp%pos (1,iparticle)),& - ABS(particle_set(iparticle)%r(2)-tmp%pos (2,iparticle)),& - ABS(particle_set(iparticle)%r(3)-tmp%pos (3,iparticle))) - tmp%max_dvel = MAX(tmp%max_dvel,ABS(particle_set(iparticle)%v(1)-tmp%vel (1,iparticle)),& - ABS(particle_set(iparticle)%v(2)-tmp%vel (2,iparticle)),& - ABS(particle_set(iparticle)%v(3)-tmp%vel (3,iparticle))) + tmp%max_vel = MAX(tmp%max_vel,ABS(tmp%vel(1,iparticle)),& + ABS(tmp%vel (2,iparticle)),ABS(tmp%vel (3,iparticle))) + tmp%max_vel_sc = MAX(tmp%max_vel_sc,& + ABS(tmp%shell_vel(1,shell_index)-tmp%core_vel(1,shell_index)),& + ABS(tmp%shell_vel(2,shell_index)-tmp%core_vel(2,shell_index)),& + ABS(tmp%shell_vel(3,shell_index)-tmp%core_vel(3,shell_index))) + tmp%max_dr = MAX(tmp%max_dr,& + ABS(particle_set(iparticle)%r(1)-tmp%pos (1,iparticle)),& + ABS(particle_set(iparticle)%r(2)-tmp%pos (2,iparticle)),& + ABS(particle_set(iparticle)%r(3)-tmp%pos (3,iparticle))) + tmp%max_dvel = MAX(tmp%max_dvel,& + ABS(particle_set(iparticle)%v(1)-tmp%vel (1,iparticle)),& + ABS(particle_set(iparticle)%v(2)-tmp%vel (2,iparticle)),& + ABS(particle_set(iparticle)%v(3)-tmp%vel (3,iparticle))) - dsc(:) = tmp%shell_pos(:,shell_index)-tmp%core_pos(:,shell_index)- & - shell_particle_set(shell_index)%r(:)+core_particle_set(shell_index)%r(:) - tmp%max_dsc = MAX(tmp%max_dsc,ABS(dsc(1)),ABS(dsc(2)),ABS(dsc(3))) - dvsc(:) = tmp%shell_vel(:,shell_index)-tmp%core_vel(:,shell_index)- & - shell_particle_set(shell_index)%v(:)+core_particle_set(shell_index)%v(:) - tmp%max_dvel_sc = MAX(tmp%max_dvel_sc,ABS(dvsc(1)),ABS(dvsc(2)),ABS(dvsc(3))) - END DO ! iparticle_local + dsc(:) = tmp%shell_pos(:,shell_index)-tmp%core_pos(:,shell_index)- & + shell_particle_set(shell_index)%r(:)+core_particle_set(shell_index)%r(:) + tmp%max_dsc = MAX(tmp%max_dsc,ABS(dsc(1)),ABS(dsc(2)),ABS(dsc(3))) + dvsc(:) = tmp%shell_vel(:,shell_index)-tmp%core_vel(:,shell_index)- & + shell_particle_set(shell_index)%v(:)+core_particle_set(shell_index)%v(:) + tmp%max_dvel_sc = MAX(tmp%max_dvel_sc,ABS(dvsc(1)),ABS(dvsc(2)),ABS(dvsc(3))) + END DO ! iparticle_local + ELSE + DO iparticle_local=1,nparticle_local + iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) + shell_index = particle_set(iparticle)%shell_index + tmp%shell_vel (:,shell_index) = & + shell_particle_set(shell_index)%v(:) * tmp%scale_v(:) * tmp%scale_v(:) +& + tmp%scale_v(:) * tmp%poly_v(:) * dms * shell_particle_set(shell_index)%f(:) + tmp%shell_pos (:,shell_index) = & + shell_particle_set(shell_index)%r(:) * tmp%scale_r(:) * tmp%scale_r(:) +& + tmp%scale_r(:) * tmp%poly_r(:) *dt * tmp%shell_vel ( :, shell_index ) + tmp%core_vel (:,shell_index) = & + core_particle_set(shell_index)%v(:) * tmp%scale_v(:) * tmp%scale_v(:) +& + tmp%scale_v(:) * tmp%poly_v(:) *dmc * core_particle_set(shell_index)%f(:) + tmp%core_pos (:,shell_index) = & + core_particle_set(shell_index)%r(:) * tmp%scale_r(:) * tmp%scale_r(:) +& + tmp%scale_r(:) * tmp%poly_r(:) *dt * tmp%core_vel ( :, shell_index ) + + tmp%vel (:,iparticle) = fac_masss*tmp%shell_vel (:,shell_index)+& + fac_massc*tmp%core_vel (:,shell_index) + tmp%pos (:,iparticle) = fac_masss*tmp%shell_pos (:,shell_index)+& + fac_massc*tmp%core_pos (:,shell_index) + tmp%max_vel = MAX(tmp%max_vel,& + ABS(tmp%vel(1,iparticle)),ABS(tmp%vel (2,iparticle)),ABS(tmp%vel (3,iparticle))) + tmp%max_vel_sc = MAX(tmp%max_vel_sc,& + ABS(tmp%shell_vel(1,shell_index)-tmp%core_vel(1,shell_index)),& + ABS(tmp%shell_vel(2,shell_index)-tmp%core_vel(2,shell_index)),& + ABS(tmp%shell_vel(3,shell_index)-tmp%core_vel(3,shell_index))) + tmp%max_dr = MAX(tmp%max_dr,& + ABS(particle_set(iparticle)%r(1)-tmp%pos (1,iparticle)),& + ABS(particle_set(iparticle)%r(2)-tmp%pos (2,iparticle)),& + ABS(particle_set(iparticle)%r(3)-tmp%pos (3,iparticle))) + tmp%max_dvel = MAX(tmp%max_dvel,& + ABS(particle_set(iparticle)%v(1)-tmp%vel (1,iparticle)),& + ABS(particle_set(iparticle)%v(2)-tmp%vel (2,iparticle)),& + ABS(particle_set(iparticle)%v(3)-tmp%vel (3,iparticle))) + dsc(:) = tmp%shell_pos (:,shell_index)-tmp%core_pos (:,shell_index)-& + shell_particle_set(shell_index)%r(:)+core_particle_set(shell_index)%r(:) + tmp%max_dsc = MAX(tmp%max_dsc,ABS(dsc(1)),ABS(dsc(2)),ABS(dsc(3))) + dvsc(:) = tmp%shell_vel(:,shell_index)-tmp%core_vel(:,shell_index)- & + shell_particle_set(shell_index)%v(:)+core_particle_set(shell_index)%v(:) + tmp%max_dvel_sc = MAX(tmp%max_dvel_sc,ABS(dvsc(1)),ABS(dvsc(2)),ABS(dvsc(3))) + END DO ! iparticle_local + END IF ELSE - DO iparticle_local=1,nparticle_local - iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) - shell_index = particle_set(iparticle)%shell_index - tmp%shell_vel (:,shell_index) = & - shell_particle_set(shell_index)%v(:) * tmp%scale_v(:) * tmp%scale_v(:) +& - tmp%scale_v(:) * tmp%poly_v(:) * dms * shell_particle_set(shell_index)%f(:) - tmp%shell_pos (:,shell_index) = & - shell_particle_set(shell_index)%r(:) * tmp%scale_r(:) * tmp%scale_r(:) +& - tmp%scale_r(:) * tmp%poly_r(:) *dt * tmp%shell_vel ( :, shell_index ) - tmp%core_vel (:,shell_index) = & - core_particle_set(shell_index)%v(:) * tmp%scale_v(:) * tmp%scale_v(:) +& - tmp%scale_v(:) * tmp%poly_v(:) *dmc * core_particle_set(shell_index)%f(:) - tmp%core_pos (:,shell_index) = & - core_particle_set(shell_index)%r(:) * tmp%scale_r(:) * tmp%scale_r(:) +& - tmp%scale_r(:) * tmp%poly_r(:) *dt * tmp%core_vel ( :, shell_index ) - - tmp%vel (:,iparticle) = fac_masss*tmp%shell_vel (:,shell_index)+& - fac_massc*tmp%core_vel (:,shell_index) - tmp%pos (:,iparticle) = fac_masss*tmp%shell_pos (:,shell_index)+& - fac_massc*tmp%core_pos (:,shell_index) - tmp%max_vel = MAX(tmp%max_vel,ABS(tmp%vel(1,iparticle)),ABS(tmp%vel (2,iparticle)),ABS(tmp%vel (3,iparticle))) - tmp%max_vel_sc = MAX(tmp%max_vel_sc,ABS(tmp%shell_vel(1,shell_index)-tmp%core_vel(1,shell_index)),& - ABS(tmp%shell_vel(2,shell_index)-tmp%core_vel(2,shell_index)),& - ABS(tmp%shell_vel(3,shell_index)-tmp%core_vel(3,shell_index))) - tmp%max_dr = MAX(tmp%max_dr,ABS(particle_set(iparticle)%r(1)-tmp%pos (1,iparticle)),& - ABS(particle_set(iparticle)%r(2)-tmp%pos (2,iparticle)),& - ABS(particle_set(iparticle)%r(3)-tmp%pos (3,iparticle))) - tmp%max_dvel = MAX(tmp%max_dvel,ABS(particle_set(iparticle)%v(1)-tmp%vel (1,iparticle)),& - ABS(particle_set(iparticle)%v(2)-tmp%vel (2,iparticle)),& - ABS(particle_set(iparticle)%v(3)-tmp%vel (3,iparticle))) - dsc(:) = tmp%shell_pos (:,shell_index)-tmp%core_pos (:,shell_index)-& - shell_particle_set(shell_index)%r(:)+core_particle_set(shell_index)%r(:) - tmp%max_dsc = MAX(tmp%max_dsc,ABS(dsc(1)),ABS(dsc(2)),ABS(dsc(3))) - dvsc(:) = tmp%shell_vel(:,shell_index)-tmp%core_vel(:,shell_index)- & - shell_particle_set(shell_index)%v(:)+core_particle_set(shell_index)%v(:) - tmp%max_dvel_sc = MAX(tmp%max_dvel_sc,ABS(dvsc(1)),ABS(dvsc(2)),ABS(dvsc(3))) - END DO ! iparticle_local - END IF - ELSE - nparticle_local = local_particles%n_el(iparticle_kind) - IF(PRESENT(u)) THEN - DO iparticle_local=1,nparticle_local - iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) - ! Transform positions and velocities and forces - CALL transform_first(particle_set,tmp%pos,tmp%vel,& - iparticle,u,dm,dt,tmp%poly_v,tmp%poly_r,tmp%scale_v,tmp%scale_r) - - tmp%max_vel = MAX(tmp%max_vel,ABS(tmp%vel (1,iparticle)),ABS(tmp%vel (2,iparticle)),ABS(tmp%vel (3,iparticle))) - tmp%max_dr = MAX(tmp%max_dr,ABS(particle_set(iparticle)%r(1)-tmp%pos (1,iparticle)),& - ABS(particle_set(iparticle)%r(2)-tmp%pos (2,iparticle)),& - ABS(particle_set(iparticle)%r(3)-tmp%pos (3,iparticle))) - tmp%max_dvel = MAX(tmp%max_dvel,ABS(particle_set(iparticle)%v(1)-tmp%vel (1,iparticle)),& - ABS(particle_set(iparticle)%v(2)-tmp%vel (2,iparticle)),& - ABS(particle_set(iparticle)%v(3)-tmp%vel (3,iparticle))) - END DO ! iparticle_local - ELSE - DO iparticle_local=1,nparticle_local - iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) - tmp%vel (1,iparticle) = & - particle_set(iparticle)%v(1) * tmp%scale_v(1) * tmp%scale_v(1)+& - tmp%scale_v(1) * tmp%poly_v(1) * dm * particle_set(iparticle)%f(1) - tmp%vel (2,iparticle) = & - particle_set(iparticle)%v(2) * tmp%scale_v(2) * tmp%scale_v(2)+& - tmp%scale_v(2) * tmp%poly_v(2) * dm * particle_set(iparticle)%f(2) - tmp%vel (3,iparticle) = & - particle_set(iparticle)%v(3) * tmp%scale_v(3) * tmp%scale_v(3)+& - tmp%scale_v(3) * tmp%poly_v(3) * dm * particle_set(iparticle)%f(3) - tmp%pos (1,iparticle) = & - particle_set(iparticle)%r(1) * tmp%scale_r(1) * tmp%scale_r(1)+& - tmp%scale_r(1) * tmp%poly_r(1) * dt * tmp%vel ( 1, iparticle ) - tmp%pos (2,iparticle) = & - particle_set(iparticle)%r(2) * tmp%scale_r(2) * tmp%scale_r(2)+& - tmp%scale_r(2) * tmp%poly_r(2) * dt * tmp%vel ( 2, iparticle ) - tmp%pos (3,iparticle) = & - particle_set(iparticle)%r(3) * tmp%scale_r(3) * tmp%scale_r(3)+& - tmp%scale_r(3) * tmp%poly_r(3) * dt * tmp%vel ( 3, iparticle ) - - tmp%max_vel = MAX(tmp%max_vel,ABS(tmp%vel (1,iparticle)),ABS(tmp%vel (2,iparticle)),ABS(tmp%vel (3,iparticle))) - tmp%max_dr = MAX(tmp%max_dr,ABS(particle_set(iparticle)%r(1)-tmp%pos (1,iparticle)),& - ABS(particle_set(iparticle)%r(2)-tmp%pos (2,iparticle)),& - ABS(particle_set(iparticle)%r(3)-tmp%pos (3,iparticle))) - tmp%max_dvel = MAX(tmp%max_dvel,ABS(particle_set(iparticle)%v(1)-tmp%vel (1,iparticle)),& - ABS(particle_set(iparticle)%v(2)-tmp%vel (2,iparticle)),& - ABS(particle_set(iparticle)%v(3)-tmp%vel (3,iparticle))) - END DO ! iparticle_local + nparticle_local = local_particles%n_el(iparticle_kind) + IF(PRESENT(u)) THEN + DO iparticle_local=1,nparticle_local + iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) + ! Transform positions and velocities and forces + CALL transform_first(particle_set,tmp%pos,tmp%vel,& + iparticle,u,dm,dt,tmp%poly_v,tmp%poly_r,tmp%scale_v,tmp%scale_r) + + tmp%max_vel = MAX(tmp%max_vel,& + ABS(tmp%vel (1,iparticle)),ABS(tmp%vel (2,iparticle)),ABS(tmp%vel (3,iparticle))) + tmp%max_dr = MAX(tmp%max_dr,ABS(particle_set(iparticle)%r(1)-tmp%pos (1,iparticle)),& + ABS(particle_set(iparticle)%r(2)-tmp%pos (2,iparticle)),& + ABS(particle_set(iparticle)%r(3)-tmp%pos (3,iparticle))) + tmp%max_dvel = MAX(tmp%max_dvel,& + ABS(particle_set(iparticle)%v(1)-tmp%vel (1,iparticle)),& + ABS(particle_set(iparticle)%v(2)-tmp%vel (2,iparticle)),& + ABS(particle_set(iparticle)%v(3)-tmp%vel (3,iparticle))) + END DO ! iparticle_local + ELSE + DO iparticle_local=1,nparticle_local + iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) + tmp%vel (1,iparticle) = & + particle_set(iparticle)%v(1) * tmp%scale_v(1) * tmp%scale_v(1)+& + tmp%scale_v(1) * tmp%poly_v(1) * dm * particle_set(iparticle)%f(1) + tmp%vel (2,iparticle) = & + particle_set(iparticle)%v(2) * tmp%scale_v(2) * tmp%scale_v(2)+& + tmp%scale_v(2) * tmp%poly_v(2) * dm * particle_set(iparticle)%f(2) + tmp%vel (3,iparticle) = & + particle_set(iparticle)%v(3) * tmp%scale_v(3) * tmp%scale_v(3)+& + tmp%scale_v(3) * tmp%poly_v(3) * dm * particle_set(iparticle)%f(3) + tmp%pos (1,iparticle) = & + particle_set(iparticle)%r(1) * tmp%scale_r(1) * tmp%scale_r(1)+& + tmp%scale_r(1) * tmp%poly_r(1) * dt * tmp%vel ( 1, iparticle ) + tmp%pos (2,iparticle) = & + particle_set(iparticle)%r(2) * tmp%scale_r(2) * tmp%scale_r(2)+& + tmp%scale_r(2) * tmp%poly_r(2) * dt * tmp%vel ( 2, iparticle ) + tmp%pos (3,iparticle) = & + particle_set(iparticle)%r(3) * tmp%scale_r(3) * tmp%scale_r(3)+& + tmp%scale_r(3) * tmp%poly_r(3) * dt * tmp%vel ( 3, iparticle ) + + tmp%max_vel = MAX(tmp%max_vel,& + ABS(tmp%vel (1,iparticle)),ABS(tmp%vel (2,iparticle)),ABS(tmp%vel (3,iparticle))) + tmp%max_dr = MAX(tmp%max_dr,& + ABS(particle_set(iparticle)%r(1)-tmp%pos (1,iparticle)),& + ABS(particle_set(iparticle)%r(2)-tmp%pos (2,iparticle)),& + ABS(particle_set(iparticle)%r(3)-tmp%pos (3,iparticle))) + tmp%max_dvel = MAX(tmp%max_dvel,& + ABS(particle_set(iparticle)%v(1)-tmp%vel (1,iparticle)),& + ABS(particle_set(iparticle)%v(2)-tmp%vel (2,iparticle)),& + ABS(particle_set(iparticle)%v(3)-tmp%vel (3,iparticle))) + END DO + END IF END IF END IF - END DO ! iparticle_kind - + END DO END SUBROUTINE vv_first ! ***************************************************************************** @@ -1191,90 +1208,91 @@ CONTAINS atomic_kind => atomic_kind_set(iparticle_kind) CALL get_atomic_kind(atomic_kind=atomic_kind,mass=mass,& shell_active=is_shell) - dm = 0.5_dp*dt/mass - IF(is_shell .AND. shell_adiabatic) THEN - CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell) - dms = 0.5_dp*dt/shell%mass_shell - dmc = 0.5_dp*dt/shell%mass_core - fac_masss = shell%mass_shell/mass - fac_massc = shell%mass_core/mass - nparticle_local = local_particles%n_el(iparticle_kind) - IF(PRESENT(u)) THEN - DO iparticle_local=1,nparticle_local - iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) - shell_index = particle_set(iparticle)%shell_index - ! Transform velocities and forces shell - CALL transform_second(shell_particle_set,tmp%shell_vel,shell_index,& - u,dms,tmp%poly_v,tmp%scale_v) + IF (mass/=0.0_dp) THEN + dm = 0.5_dp*dt/mass + IF(is_shell .AND. shell_adiabatic) THEN + CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell) + dms = 0.5_dp*dt/shell%mass_shell + dmc = 0.5_dp*dt/shell%mass_core + fac_masss = shell%mass_shell/mass + fac_massc = shell%mass_core/mass + nparticle_local = local_particles%n_el(iparticle_kind) + IF(PRESENT(u)) THEN + DO iparticle_local=1,nparticle_local + iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) + shell_index = particle_set(iparticle)%shell_index + ! Transform velocities and forces shell + CALL transform_second(shell_particle_set,tmp%shell_vel,shell_index,& + u,dms,tmp%poly_v,tmp%scale_v) - ! Transform velocities and forces core - CALL transform_second(core_particle_set,tmp%core_vel,shell_index,& - u,dmc,tmp%poly_v,tmp%scale_v) + ! Transform velocities and forces core + CALL transform_second(core_particle_set,tmp%core_vel,shell_index,& + u,dmc,tmp%poly_v,tmp%scale_v) - ! Derive velocties of the COM - tmp%vel(1,iparticle) = fac_masss*tmp%shell_vel(1,shell_index) + & - fac_massc*tmp%core_vel(1,shell_index) - tmp%vel(2,iparticle) = fac_masss*tmp%shell_vel(2,shell_index)+ & - fac_massc*tmp%core_vel(2,shell_index) - tmp%vel(3,iparticle) = fac_masss*tmp%shell_vel(3,shell_index) + & - fac_massc*tmp%core_vel(3,shell_index) - END DO ! iparticle_local + ! Derive velocties of the COM + tmp%vel(1,iparticle) = fac_masss*tmp%shell_vel(1,shell_index) + & + fac_massc*tmp%core_vel(1,shell_index) + tmp%vel(2,iparticle) = fac_masss*tmp%shell_vel(2,shell_index)+ & + fac_massc*tmp%core_vel(2,shell_index) + tmp%vel(3,iparticle) = fac_masss*tmp%shell_vel(3,shell_index) + & + fac_massc*tmp%core_vel(3,shell_index) + END DO ! iparticle_local + ELSE + DO iparticle_local=1,nparticle_local + iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) + shell_index = particle_set(iparticle)%shell_index + tmp%shell_vel (1,shell_index) = & + tmp%shell_vel (1,shell_index) * tmp%scale_v(1) * tmp%scale_v(1) + & + tmp%scale_v(1) * tmp%poly_v(1) * dms*shell_particle_set(shell_index)%f(1) + tmp%shell_vel (2,shell_index) =& + tmp%shell_vel (2,shell_index) * tmp%scale_v(2) * tmp%scale_v(2) + & + tmp%scale_v(2) * tmp%poly_v(2) * dms*shell_particle_set(shell_index)%f(2) + tmp%shell_vel (3,shell_index) =& + tmp%shell_vel (3,shell_index) * tmp%scale_v(3) * tmp%scale_v(3) + & + tmp%scale_v(3) * tmp%poly_v(3) * dms*shell_particle_set(shell_index)%f(3) + tmp%core_vel (1,shell_index) =& + tmp%core_vel (1,shell_index) * tmp%scale_v(1) * tmp%scale_v(1) + & + tmp%scale_v(1) * tmp%poly_v(1) * dmc*core_particle_set(shell_index)%f(1) + tmp%core_vel (2,shell_index) =& + tmp%core_vel (2,shell_index) * tmp%scale_v(2) * tmp%scale_v(2)+ & + tmp%scale_v(2) * tmp%poly_v(2) * dmc*core_particle_set(shell_index)%f(2) + tmp%core_vel (3,shell_index) =& + tmp%core_vel (3,shell_index) * tmp%scale_v(3) * tmp%scale_v(3)+ & + tmp%scale_v(3) * tmp%poly_v(3) * dmc*core_particle_set(shell_index)%f(3) + + tmp%vel ( 1, iparticle) = fac_masss*tmp%shell_vel (1,shell_index)+ & + fac_massc*tmp%core_vel (1,shell_index) + tmp%vel ( 2, iparticle) = fac_masss*tmp%shell_vel (2,shell_index)+ & + fac_massc*tmp%core_vel (2,shell_index) + tmp%vel ( 3, iparticle) = fac_masss*tmp%shell_vel (3,shell_index)+ & + fac_massc*tmp%core_vel (3,shell_index) + END DO ! iparticle_local + END IF ELSE - DO iparticle_local=1,nparticle_local - iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) - shell_index = particle_set(iparticle)%shell_index - tmp%shell_vel (1,shell_index) = & - tmp%shell_vel (1,shell_index) * tmp%scale_v(1) * tmp%scale_v(1) + & - tmp%scale_v(1) * tmp%poly_v(1) * dms*shell_particle_set(shell_index)%f(1) - tmp%shell_vel (2,shell_index) =& - tmp%shell_vel (2,shell_index) * tmp%scale_v(2) * tmp%scale_v(2) + & - tmp%scale_v(2) * tmp%poly_v(2) * dms*shell_particle_set(shell_index)%f(2) - tmp%shell_vel (3,shell_index) =& - tmp%shell_vel (3,shell_index) * tmp%scale_v(3) * tmp%scale_v(3) + & - tmp%scale_v(3) * tmp%poly_v(3) * dms*shell_particle_set(shell_index)%f(3) - tmp%core_vel (1,shell_index) =& - tmp%core_vel (1,shell_index) * tmp%scale_v(1) * tmp%scale_v(1) + & - tmp%scale_v(1) * tmp%poly_v(1) * dmc*core_particle_set(shell_index)%f(1) - tmp%core_vel (2,shell_index) =& - tmp%core_vel (2,shell_index) * tmp%scale_v(2) * tmp%scale_v(2)+ & - tmp%scale_v(2) * tmp%poly_v(2) * dmc*core_particle_set(shell_index)%f(2) - tmp%core_vel (3,shell_index) =& - tmp%core_vel (3,shell_index) * tmp%scale_v(3) * tmp%scale_v(3)+ & - tmp%scale_v(3) * tmp%poly_v(3) * dmc*core_particle_set(shell_index)%f(3) - - tmp%vel ( 1, iparticle) = fac_masss*tmp%shell_vel (1,shell_index)+ & - fac_massc*tmp%core_vel (1,shell_index) - tmp%vel ( 2, iparticle) = fac_masss*tmp%shell_vel (2,shell_index)+ & - fac_massc*tmp%core_vel (2,shell_index) - tmp%vel ( 3, iparticle) = fac_masss*tmp%shell_vel (3,shell_index)+ & - fac_massc*tmp%core_vel (3,shell_index) - - END DO ! iparticle_local - END IF - ELSE - nparticle_local = local_particles%n_el(iparticle_kind) - IF(PRESENT(u)) THEN - DO iparticle_local=1,nparticle_local - iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) - CALL transform_second(particle_set,tmp%vel,iparticle,& - u,dm,tmp%poly_v,tmp%scale_v) - END DO ! iparticle_local - ELSE - DO iparticle_local=1,nparticle_local - iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) - tmp%vel ( 1, iparticle) = & - tmp%vel ( 1, iparticle) * tmp%scale_v(1) * tmp%scale_v(1) + & - tmp%scale_v(1) * tmp%poly_v(1) * dm*particle_set(iparticle)%f(1) - tmp%vel ( 2, iparticle) = & - tmp%vel ( 2, iparticle) * tmp%scale_v(2) * tmp%scale_v(2) + & - tmp%scale_v(2) * tmp%poly_v(2) * dm*particle_set(iparticle)%f(2) - tmp%vel ( 3, iparticle) = & - tmp%vel ( 3, iparticle) * tmp%scale_v(3) * tmp%scale_v(3) + & - tmp%scale_v(3) * tmp%poly_v(3) * dm*particle_set(iparticle)%f(3) - END DO ! iparticle_local + nparticle_local = local_particles%n_el(iparticle_kind) + IF(PRESENT(u)) THEN + DO iparticle_local=1,nparticle_local + iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) + CALL transform_second(particle_set,tmp%vel,iparticle,& + u,dm,tmp%poly_v,tmp%scale_v) + END DO ! iparticle_local + ELSE + DO iparticle_local=1,nparticle_local + iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) + tmp%vel ( 1, iparticle) = & + tmp%vel ( 1, iparticle) * tmp%scale_v(1) * tmp%scale_v(1) + & + tmp%scale_v(1) * tmp%poly_v(1) * dm*particle_set(iparticle)%f(1) + tmp%vel ( 2, iparticle) = & + tmp%vel ( 2, iparticle) * tmp%scale_v(2) * tmp%scale_v(2) + & + tmp%scale_v(2) * tmp%poly_v(2) * dm*particle_set(iparticle)%f(2) + tmp%vel ( 3, iparticle) = & + tmp%vel ( 3, iparticle) * tmp%scale_v(3) * tmp%scale_v(3) + & + tmp%scale_v(3) * tmp%poly_v(3) * dm*particle_set(iparticle)%f(3) + END DO + END IF END IF END IF - END DO ! iparticle_kind + END DO END SUBROUTINE vv_second diff --git a/src/md_vel_utils.F b/src/md_vel_utils.F index 23e78e5c31..f341273a1a 100644 --- a/src/md_vel_utils.F +++ b/src/md_vel_utils.F @@ -9,6 +9,8 @@ !> \author CJM !> \author Teodoro Laino [tlaino] - University of Zurich - 10.2008 !> reorganization of the original routines/modules +!> Teodoro Laino [tlaino] 12.2008 - Preparing for VIRTUAL SITE constraints +!> (patch by Marcel Baer) ! ***************************************************************************** MODULE md_vel_utils USE atomic_kind_list_types, ONLY: atomic_kind_list_type @@ -142,13 +144,15 @@ CONTAINS DO i = 1, SIZE(part) atomic_kind => part(i)%atomic_kind CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) - SELECT CASE(is_fixed(i)) - CASE(use_perd_x,use_perd_y,use_perd_z,use_perd_xy,use_perd_xz,use_perd_yz,use_perd_none) - vcom(1) = vcom(1) + part(i)%v(1) * mass - vcom(2) = vcom(2) + part(i)%v(2) * mass - vcom(3) = vcom(3) + part(i)%v(3) * mass - denom = denom + mass - END SELECT + IF (mass.NE.0.0) THEN + SELECT CASE(is_fixed(i)) + CASE(use_perd_x,use_perd_y,use_perd_z,use_perd_xy,use_perd_xz,use_perd_yz,use_perd_none) + vcom(1) = vcom(1) + part(i)%v(1) * mass + vcom(2) = vcom(2) + part(i)%v(2) * mass + vcom(3) = vcom(3) + part(i)%v(3) * mass + denom = denom + mass + END SELECT + END IF END DO vcom = vcom/denom IF (PRESENT(ecom)) THEN @@ -707,27 +711,29 @@ CONTAINS part(i)%v(1) = 0.0_dp part(i)%v(2) = 0.0_dp part(i)%v(3) = 0.0_dp - SELECT CASE(is_fixed(i)) - CASE (use_perd_x) - part(i)%v(2) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) - part(i)%v(3) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) - CASE (use_perd_y) - part(i)%v(1) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) - part(i)%v(3) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) - CASE (use_perd_z) - part(i)%v(1) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) - part(i)%v(2) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) - CASE (use_perd_xy) - part(i)%v(3) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) - CASE (use_perd_xz) - part(i)%v(2) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) - CASE (use_perd_yz) - part(i)%v(1) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) - CASE (use_perd_none) - part(i)%v(1) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) - part(i)%v(2) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) - part(i)%v(3) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) - END SELECT + IF (mass.NE.0.0) THEN + SELECT CASE(is_fixed(i)) + CASE (use_perd_x) + part(i)%v(2) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) + part(i)%v(3) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) + CASE (use_perd_y) + part(i)%v(1) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) + part(i)%v(3) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) + CASE (use_perd_z) + part(i)%v(1) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) + part(i)%v(2) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) + CASE (use_perd_xy) + part(i)%v(3) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) + CASE (use_perd_xz) + part(i)%v(2) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) + CASE (use_perd_yz) + part(i)%v(1) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) + CASE (use_perd_none) + part(i)%v(1) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) + part(i)%v(2) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) + part(i)%v(3) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) + END SELECT + END IF END DO ! Subtract the vcom diff --git a/src/molecule_kind_types.F b/src/molecule_kind_types.F index 5c1efdca8e..8b9228c187 100644 --- a/src/molecule_kind_types.F +++ b/src/molecule_kind_types.F @@ -7,7 +7,8 @@ !> \brief Define the molecule kind structure types and the corresponding !> functionality !> \par History -!> none +!> Teodoro Laino [tlaino] 12.2008 - Preparing for VIRTUAL SITE constraints +!> (patch by Marcel Baer) !> \author MK (22.08.2003) ! ***************************************************************************** MODULE molecule_kind_types @@ -138,6 +139,13 @@ MODULE molecule_kind_types TYPE(restraint_type) :: restraint END TYPE g4x6_constraint_type +! ***************************************************************************** + TYPE vsite_constraint_type + INTEGER :: a,b,c,d + REAL(KIND = dp) :: wbc, wdc + TYPE(restraint_type) :: restraint + END TYPE vsite_constraint_type + ! ***************************************************************************** TYPE fixd_constraint_type TYPE(restraint_type) :: restraint @@ -166,9 +174,10 @@ MODULE molecule_kind_types TYPE(onfo_kind_type), DIMENSION(:), POINTER :: onfo_kind_set TYPE(onfo_type), DIMENSION(:), POINTER :: onfo_list TYPE(colvar_constraint_type), DIMENSION(:), POINTER :: colv_list - TYPE(g3x3_constraint_type), DIMENSION(:), POINTER :: g3x3_list - TYPE(g4x6_constraint_type), DIMENSION(:), POINTER :: g4x6_list - TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list + TYPE(g3x3_constraint_type), DIMENSION(:), POINTER :: g3x3_list + TYPE(g4x6_constraint_type), DIMENSION(:), POINTER :: g4x6_list + TYPE(vsite_constraint_type), DIMENSION(:), POINTER :: vsite_list + TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list TYPE(shell_type), DIMENSION(:), POINTER :: shell_list CHARACTER(LEN=default_string_length) :: name REAL(KIND = dp) :: charge,& @@ -183,6 +192,7 @@ MODULE molecule_kind_types nub,& ng3x3, ng3x3_restraint,& ng4x6, ng4x6_restraint,& + nvsite, nvsite_restraint,& nfixd, nfixd_restraint,& nmolecule, nshell TYPE(colvar_counters) :: ncolv @@ -214,6 +224,7 @@ MODULE molecule_kind_types colvar_constraint_type,& g3x3_constraint_type,& g4x6_constraint_type,& + vsite_constraint_type,& fixd_constraint_type,& local_fixd_constraint_type,& molecule_kind_type,& @@ -327,6 +338,7 @@ CONTAINS NULLIFY (molecule_kind_set(imolecule_kind)%onfo_kind_set) NULLIFY (molecule_kind_set(imolecule_kind)%g3x3_list) NULLIFY (molecule_kind_set(imolecule_kind)%g4x6_list) + NULLIFY (molecule_kind_set(imolecule_kind)%vsite_list) NULLIFY (molecule_kind_set(imolecule_kind)%fixd_list) NULLIFY (molecule_kind_set(imolecule_kind)%shell_list) NULLIFY (molecule_kind_set(imolecule_kind)%torsion_list) @@ -348,9 +360,11 @@ CONTAINS molecule_kind_set(imolecule_kind)%ncolv,error) molecule_kind_set(imolecule_kind)%ng3x3 = 0 molecule_kind_set(imolecule_kind)%ng4x6 = 0 + molecule_kind_set(imolecule_kind)%nvsite = 0 molecule_kind_set(imolecule_kind)%nfixd = 0 molecule_kind_set(imolecule_kind)%ng3x3_restraint = 0 molecule_kind_set(imolecule_kind)%ng4x6_restraint = 0 + molecule_kind_set(imolecule_kind)%nvsite_restraint = 0 molecule_kind_set(imolecule_kind)%nfixd_restraint = 0 molecule_kind_set(imolecule_kind)%nmolecule = 0 molecule_kind_set(imolecule_kind)%ntorsion = 0 @@ -447,6 +461,10 @@ CONTAINS DEALLOCATE (molecule_kind_set(imolecule_kind)%g4x6_list,STAT=stat) CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) END IF + IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%vsite_list)) THEN + DEALLOCATE (molecule_kind_set(imolecule_kind)%vsite_list,STAT=stat) + CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) + END IF IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%fixd_list)) THEN DEALLOCATE (molecule_kind_set(imolecule_kind)%fixd_list,STAT=stat) CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) @@ -489,11 +507,11 @@ CONTAINS ! ***************************************************************************** SUBROUTINE get_molecule_kind(molecule_kind,atom_list,bond_list,bend_list,& ub_list,impr_list,onfo_list,colv_list,fixd_list,& - g3x3_list,g4x6_list,torsion_list,shell_list,& + g3x3_list,g4x6_list,vsite_list,torsion_list,shell_list,& name,mass,charge,kind_number,natom,nbend,nbond,nub,& nimpr,nonfo,nconstraint,nconstraint_fixd,nfixd,ncolv,ng3x3,ng4x6,& - nfixd_restraint,ng3x3_restraint, ng4x6_restraint,& - nrestraints,nmolecule,nsgf,nshell,ntorsion,& + nvsite,nfixd_restraint,ng3x3_restraint, ng4x6_restraint,& + nvsite_restraint,nrestraints,nmolecule,nsgf,nshell,ntorsion,& molecule_list,nelectron, bond_kind_set,bend_kind_set,& ub_kind_set,impr_kind_set,onfo_kind_set,torsion_kind_set,& molname_generated) @@ -519,6 +537,8 @@ CONTAINS DIMENSION(:), OPTIONAL, POINTER :: g3x3_list TYPE(g4x6_constraint_type), & DIMENSION(:), OPTIONAL, POINTER :: g4x6_list + TYPE(vsite_constraint_type), & + DIMENSION(:), OPTIONAL, POINTER :: vsite_list TYPE(torsion_type), DIMENSION(:), & OPTIONAL, POINTER :: torsion_list TYPE(shell_type), DIMENSION(:), & @@ -532,9 +552,9 @@ CONTAINS nconstraint_fixd, nfixd TYPE(colvar_counters), INTENT(out), & OPTIONAL :: ncolv - INTEGER, INTENT(OUT), OPTIONAL :: ng3x3, ng4x6, nfixd_restraint, & - ng3x3_restraint, ng4x6_restraint, nrestraints, nmolecule, nsgf, nshell, & - ntorsion + INTEGER, INTENT(OUT), OPTIONAL :: ng3x3, ng4x6, nvsite, nfixd_restraint, & + ng3x3_restraint, ng4x6_restraint, nvsite_restraint, nrestraints, & + nmolecule, nsgf, nshell, ntorsion INTEGER, DIMENSION(:), OPTIONAL, POINTER :: molecule_list INTEGER, INTENT(OUT), OPTIONAL :: nelectron TYPE(bond_kind_type), DIMENSION(:), & @@ -573,6 +593,7 @@ CONTAINS IF (PRESENT(colv_list)) colv_list => molecule_kind%colv_list IF (PRESENT(g3x3_list)) g3x3_list => molecule_kind%g3x3_list IF (PRESENT(g4x6_list)) g4x6_list => molecule_kind%g4x6_list + IF (PRESENT(vsite_list)) vsite_list => molecule_kind%vsite_list IF (PRESENT(fixd_list)) fixd_list => molecule_kind%fixd_list IF (PRESENT(torsion_list)) torsion_list => molecule_kind%torsion_list IF (PRESENT(shell_list)) shell_list => molecule_kind%shell_list @@ -589,10 +610,12 @@ CONTAINS IF (PRESENT(nonfo)) nonfo = molecule_kind%nonfo IF (PRESENT(nconstraint)) nconstraint = (molecule_kind%ncolv%ntot - molecule_kind%ncolv%nrestraint) +& 3*(molecule_kind%ng3x3-molecule_kind%ng3x3_restraint) +& - 6*(molecule_kind%ng4x6-molecule_kind%ng4x6_restraint) + 6*(molecule_kind%ng4x6-molecule_kind%ng4x6_restraint)+& + 3*(molecule_kind%nvsite-molecule_kind%nvsite_restraint) IF (PRESENT(ncolv)) ncolv = molecule_kind%ncolv IF (PRESENT(ng3x3)) ng3x3= molecule_kind%ng3x3 IF (PRESENT(ng4x6)) ng4x6= molecule_kind%ng4x6 + IF (PRESENT(nvsite)) nvsite= molecule_kind%nvsite ! Number of atoms that have one or more components fixed IF (PRESENT(nfixd)) nfixd= molecule_kind%nfixd ! Number of degrees of freedom fixed @@ -614,10 +637,12 @@ CONTAINS END IF IF (PRESENT(ng3x3_restraint)) ng3x3_restraint= molecule_kind%ng3x3_restraint IF (PRESENT(ng4x6_restraint)) ng4x6_restraint= molecule_kind%ng4x6_restraint + IF (PRESENT(nvsite_restraint)) nvsite_restraint= molecule_kind%nvsite_restraint IF (PRESENT(nfixd_restraint)) nfixd_restraint= molecule_kind%nfixd_restraint IF (PRESENT(nrestraints)) nrestraints = molecule_kind%ncolv%nrestraint +& molecule_kind%ng3x3_restraint +& - molecule_kind%ng4x6_restraint + molecule_kind%ng4x6_restraint+& + molecule_kind%nvsite_restraint IF (PRESENT(nmolecule)) nmolecule = molecule_kind%nmolecule IF (PRESENT(nshell)) nshell= molecule_kind%nshell IF (PRESENT(ntorsion)) ntorsion= molecule_kind%ntorsion @@ -802,8 +827,8 @@ CONTAINS nbend,bend_list,nub,ub_list,nimpr,impr_list,& nonfo,onfo_list,ntorsion,torsion_list,fixd_list,& ncolv,colv_list,ng3x3,g3x3_list,ng4x6,nfixd,g4x6_list,& - ng3x3_restraint,ng4x6_restraint,nfixd_restraint,& - nshell, shell_list,& + nvsite,vsite_list,ng3x3_restraint,ng4x6_restraint,nfixd_restraint,& + nshell, shell_list,nvsite_restraint,& bond_kind_set,bend_kind_set,ub_kind_set,& torsion_kind_set,impr_kind_set,onfo_kind_set,& nelectron,nsgf,molname_generated) @@ -845,11 +870,15 @@ CONTAINS INTEGER, INTENT(IN), OPTIONAL :: ng4x6, nfixd TYPE(g4x6_constraint_type), & DIMENSION(:), OPTIONAL, POINTER :: g4x6_list + INTEGER, INTENT(IN), OPTIONAL :: nvsite + TYPE(vsite_constraint_type), & + DIMENSION(:), OPTIONAL, POINTER :: vsite_list INTEGER, INTENT(IN), OPTIONAL :: ng3x3_restraint, & ng4x6_restraint, & nfixd_restraint, nshell TYPE(shell_type), DIMENSION(:), & OPTIONAL, POINTER :: shell_list + INTEGER, INTENT(IN), OPTIONAL :: nvsite_restraint TYPE(bond_kind_type), DIMENSION(:), & OPTIONAL, POINTER :: bond_kind_set TYPE(bend_kind_type), DIMENSION(:), & @@ -901,11 +930,14 @@ CONTAINS IF (PRESENT(ng3x3)) molecule_kind%ng3x3 = ng3x3 IF (PRESENT(g3x3_list)) molecule_kind%g3x3_list => g3x3_list IF (PRESENT(ng4x6)) molecule_kind%ng4x6 = ng4x6 + IF (PRESENT(nvsite)) molecule_kind%nvsite = nvsite IF (PRESENT(nfixd)) molecule_kind%nfixd = nfixd IF (PRESENT(nfixd_restraint)) molecule_kind%nfixd_restraint = nfixd_restraint IF (PRESENT(ng3x3_restraint)) molecule_kind%ng3x3_restraint = ng3x3_restraint IF (PRESENT(ng4x6_restraint)) molecule_kind%ng4x6_restraint = ng4x6_restraint + IF (PRESENT(nvsite_restraint)) molecule_kind%nvsite_restraint = nvsite_restraint IF (PRESENT(g4x6_list)) molecule_kind%g4x6_list => g4x6_list + IF (PRESENT(vsite_list)) molecule_kind%vsite_list => vsite_list IF (PRESENT(fixd_list)) molecule_kind%fixd_list => fixd_list IF (PRESENT(bond_kind_set)) molecule_kind%bond_kind_set => bond_kind_set IF (PRESENT(bend_kind_set)) molecule_kind%bend_kind_set => bend_kind_set diff --git a/src/molecule_types_new.F b/src/molecule_types_new.F index af2553d2a2..ab66d01883 100644 --- a/src/molecule_types_new.F +++ b/src/molecule_types_new.F @@ -7,6 +7,8 @@ !> \brief Define the data structure for the molecule information. !> \par History !> 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) ! ***************************************************************************** MODULE molecule_types_new @@ -21,7 +23,8 @@ MODULE molecule_types_new fixd_constraint_type,& g3x3_constraint_type,& g4x6_constraint_type,& - molecule_kind_type + molecule_kind_type,& + vsite_constraint_type USE termination, ONLY: stop_memory,& stop_program #include "cp_common_uses.h" @@ -86,10 +89,12 @@ MODULE molecule_types_new INTEGER :: ntot, nrestraint INTEGER :: ng3x3, ng3x3_restraint INTEGER :: ng4x6, ng4x6_restraint + INTEGER :: nvsite, nvsite_restraint TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list TYPE(colvar_constraint_type), DIMENSION(:), POINTER :: colv_list TYPE(g3x3_constraint_type), DIMENSION(:), POINTER :: g3x3_list TYPE(g4x6_constraint_type), DIMENSION(:), POINTER :: g4x6_list + TYPE(vsite_constraint_type), DIMENSION(:), POINTER :: vsite_list TYPE(local_colvar_constraint_type), DIMENSION ( : ) , POINTER :: lcolv TYPE(local_g3x3_constraint_type), DIMENSION ( : ) , POINTER :: lg3x3 TYPE(local_g4x6_constraint_type), DIMENSION ( : ) , POINTER :: lg4x6 diff --git a/src/topology_constraint_util.F b/src/topology_constraint_util.F index 2d889efd6b..bf48e1a7a8 100644 --- a/src/topology_constraint_util.F +++ b/src/topology_constraint_util.F @@ -35,7 +35,8 @@ MODULE topology_constraint_util USE molecule_kind_types, ONLY: & atom_type, bond_type, colvar_constraint_type, fixd_constraint_type, & g3x3_constraint_type, g4x6_constraint_type, get_molecule_kind, & - molecule_kind_type, set_molecule_kind, setup_colvar_counters + molecule_kind_type, set_molecule_kind, setup_colvar_counters, & + vsite_constraint_type USE molecule_types_new, ONLY: get_molecule,& global_constraint_type,& local_colvar_constraint_type,& @@ -97,7 +98,7 @@ CONTAINS k1loc, k2loc, kk, last, last_atom, m, n_start_colv, natom, nbond, & ncolv_glob, ncolv_mol, nfixd_list_gci, nfixd_restart, nfixd_restraint, & nfixed_atoms, ng3x3, ng3x3_restraint, ng4x6, ng4x6_restraint, nhdist, & - nmolecule, nrep, offset, output_unit, stat + nmolecule, nrep, nvsite, nvsite_restraint, offset, output_unit, stat INTEGER, DIMENSION(:), POINTER :: constr_x_glob, inds, & molecule_list LOGICAL :: exclude_mm, exclude_qm, failure, fix_atom_mm, & @@ -134,6 +135,8 @@ CONTAINS TYPE(molecule_type), POINTER :: molecule TYPE(section_vals_type), POINTER :: colvar_rest, fixd_restr_rest, & hbonds_section + TYPE(vsite_constraint_type), & + DIMENSION(:), POINTER :: vsite_list failure = .FALSE. NULLIFY(logger, constr_x_mol, constr_x_glob) @@ -179,15 +182,19 @@ CONTAINS gci%fixd_list,& gci%colv_list,& gci%g3x3_list,& - gci%g4x6_list) + gci%g4x6_list,& + gci%vsite_list) gci%ntot = 0 gci%ng3x3 = 0 gci%ng4x6 = 0 + gci%nvsite = 0 gci%ng3x3_restraint = 0 gci%ng4x6_restraint = 0 + gci%nvsite_restraint = 0 CALL setup_colvar_counters(gci%colv_list,gci%ncolv,error) gci%nrestraint = gci%ng3x3_restraint +& gci%ng4x6_restraint +& + gci%nvsite_restraint +& gci%ncolv%nrestraint !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- @@ -574,7 +581,54 @@ CONTAINS !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- - ! 6. Set the group fixed_atom constraint fixd_list + ! 6. Set the group vsite constraint vsite_list + !----------------------------------------------------------------------------- + IF(topology%const_vsite) THEN + CALL give_constraint_array(cons_info%const_vsite_mol,& + cons_info%const_vsite_molname,& + cons_info%vsite_intermolecular,& + constr_x_mol,& + constr_x_glob,& + molecule_kind_set,& + cons_info%vsite_exclude_qm,& + cons_info%vsite_exclude_mm,& + error) + ! Intramolecular constraints + DO ii = 1, SIZE(molecule_kind_set) + molecule_kind => molecule_kind_set(ii) + CALL get_molecule_kind(molecule_kind=molecule_kind,& + nmolecule=nmolecule,molecule_list=molecule_list) + nvsite = SIZE(constr_x_mol(ii)%constr) + ALLOCATE(vsite_list(nvsite),STAT=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + CALL setup_vsite_list(vsite_list, constr_x_mol(ii)%constr, cons_info, nvsite_restraint, error) + CALL set_molecule_kind(molecule_kind,nvsite=nvsite,nvsite_restraint=nvsite_restraint,& + vsite_list=vsite_list) + END DO + DO j = 1, SIZE(constr_x_mol) + DEALLOCATE(constr_x_mol(j)%constr,stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + END DO + DEALLOCATE(constr_x_mol,stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + ! Intermolecular constraints + IF (ASSOCIATED(constr_x_glob)) THEN + nvsite = SIZE(constr_x_glob) + ALLOCATE(vsite_list(nvsite),STAT=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + CALL setup_vsite_list(vsite_list, constr_x_glob, cons_info, nvsite_restraint, error) + gci%vsite_list => vsite_list + gci%nvsite = nvsite + gci%nvsite_restraint = nvsite_restraint + ! Total number of Intermolecular constraints + gci%ntot = gci%nvsite + gci%ntot + DEALLOCATE(constr_x_glob,stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + END IF + END IF + !----------------------------------------------------------------------------- + !----------------------------------------------------------------------------- + ! 7. Set the group fixed_atom constraint fixd_list !----------------------------------------------------------------------------- IF(topology%const_atom) THEN ALLOCATE(fixd_list_gci(SIZE(particle_set)),stat=stat) @@ -846,6 +900,7 @@ CONTAINS ! Final setup of the number of possible restraints gci%nrestraint = gci%ng3x3_restraint +& gci%ng4x6_restraint +& + gci%nvsite_restraint +& gci%ncolv%nrestraint CALL cp_print_key_finished_output(iw,logger,subsys_section,& "PRINT%TOPOLOGY_INFO/UTIL_INFO",error=error) @@ -1013,6 +1068,43 @@ CONTAINS END SUBROUTINE setup_g4x6_list ! ***************************************************************************** +!> \brief Setup the vsite_list for the packing of constraints +!> \par History +!> \author Marcel Baer [2008] +! ***************************************************************************** + SUBROUTINE setup_vsite_list(vsite_list, ilist, cons_info, nvsite_restraint, error) + TYPE(vsite_constraint_type), & + DIMENSION(:), POINTER :: vsite_list + INTEGER, DIMENSION(:), POINTER :: ilist + TYPE(constraint_info_type), POINTER :: cons_info + INTEGER, INTENT(OUT) :: nvsite_restraint + TYPE(cp_error_type), INTENT(inout) :: error + + CHARACTER(len=*), PARAMETER :: routineN = 'setup_vsite_list', & + routineP = moduleN//':'//routineN + + INTEGER :: j, nvsite + LOGICAL :: failure + + failure = .FALSE. + nvsite = 0 + nvsite_restraint = 0 + DO nvsite = 1, SIZE(ilist) + j = ilist(nvsite) + vsite_list(nvsite)%a = cons_info%const_vsite_a(j) + vsite_list(nvsite)%b = cons_info%const_vsite_b(j) + vsite_list(nvsite)%c = cons_info%const_vsite_c(j) + vsite_list(nvsite)%d = cons_info%const_vsite_d(j) + vsite_list(nvsite)%wbc = cons_info%const_vsite_wbc(j) + vsite_list(nvsite)%wdc = cons_info%const_vsite_wdc(j) + ! Restraint + vsite_list(nvsite)%restraint%active = cons_info%vsite_restraint(j) + vsite_list(nvsite)%restraint%k0 = cons_info%vsite_k0(j) + IF (vsite_list(nvsite)%restraint%active) nvsite_restraint = nvsite_restraint+ 1 + END DO + + END SUBROUTINE setup_vsite_list +! ***************************************************************************** !> \brief Setup the lcolv for the packing of constraints !> \par History !> Updated 2007 for intermolecular constraints diff --git a/src/topology_input.F b/src/topology_input.F index b31d65f1b9..2dfeeeb6ae 100644 --- a/src/topology_input.F +++ b/src/topology_input.F @@ -126,10 +126,8 @@ SUBROUTINE read_constraints_section(topology, colvar_p, constraint_section, erro LOGICAL :: explicit, failure REAL(KIND=dp), DIMENSION(:), POINTER :: rlist TYPE(constraint_info_type), POINTER :: cons_info - TYPE(section_vals_type), POINTER :: collective_section, & - fix_atom_section, & - g3x3_section, g4x6_section, & - hbonds_section + TYPE(section_vals_type), POINTER :: collective_section, fix_atom_section, & + g3x3_section, g4x6_section, hbonds_section, vsite_section cons_info => topology%cons_info failure = .FALSE. @@ -137,6 +135,7 @@ SUBROUTINE read_constraints_section(topology, colvar_p, constraint_section, erro hbonds_section => section_vals_get_subs_vals(constraint_section,"HBONDS",error=error) g3x3_section => section_vals_get_subs_vals(constraint_section,"G3X3",error=error) g4x6_section => section_vals_get_subs_vals(constraint_section,"G4X6",error=error) + vsite_section => section_vals_get_subs_vals(constraint_section,"VIRTUAL_SITE",error=error) fix_atom_section => section_vals_get_subs_vals(constraint_section,"FIXED_ATOMS",error=error) collective_section => section_vals_get_subs_vals(constraint_section,"COLLECTIVE",error=error) ! HBONDS @@ -333,6 +332,93 @@ SUBROUTINE read_constraints_section(topology, colvar_p, constraint_section, erro cons_info%const_g46_dcd(ig) = rlist(6) END DO END IF + ! virtual + CALL section_vals_get(vsite_section,explicit=explicit, n_repetition=ncons, error=error) + IF (explicit) THEN + topology%const_vsite =.TRUE. + cons_info%nconst_vsite = ncons + ! + ALLOCATE(cons_info%const_vsite_mol(ncons),stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + ALLOCATE(cons_info%const_vsite_molname(ncons),stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + ALLOCATE(cons_info%const_vsite_a(ncons),stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + ALLOCATE(cons_info%const_vsite_b(ncons),stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + ALLOCATE(cons_info%const_vsite_c(ncons),stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + ALLOCATE(cons_info%const_vsite_d(ncons),stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + ALLOCATE(cons_info%const_vsite_wbc(ncons),stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + ALLOCATE(cons_info%const_vsite_wdc(ncons),stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + ALLOCATE(cons_info%vsite_intermolecular(ncons),stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + ALLOCATE(cons_info%vsite_restraint(ncons),stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + ALLOCATE(cons_info%vsite_k0(ncons),stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + ALLOCATE(cons_info%vsite_exclude_qm(ncons),stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + ALLOCATE(cons_info%vsite_exclude_mm(ncons),stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + DO ig = 1, ncons + CALL check_restraint(vsite_section,& + is_restraint=cons_info%vsite_restraint(ig), & + k0=cons_info%vsite_k0(ig),& + i_rep_section=ig,& + label="Virtual_SITE",& + error=error) + cons_info%const_vsite_mol(ig) = 0 + cons_info%const_vsite_molname(ig) = "UNDEF" + ! Exclude QM or MM + CALL section_vals_val_get(vsite_section,"EXCLUDE_QM",i_rep_section=ig,& + l_val=cons_info%vsite_exclude_qm(ig),error=error) + CALL section_vals_val_get(vsite_section,"EXCLUDE_MM",i_rep_section=ig,& + l_val=cons_info%vsite_exclude_mm(ig),error=error) + ! Intramolecular restraint + CALL section_vals_val_get(vsite_section,"INTERMOLECULAR",i_rep_section=ig,& + l_val=cons_info%vsite_intermolecular(ig),error=error) + ! If it is intramolecular let's unset (in case user did it) + ! the molecule and molname field + IF (cons_info%vsite_intermolecular(ig)) THEN + CALL section_vals_val_unset(vsite_section,"MOLECULE",i_rep_section=ig,error=error) + CALL section_vals_val_unset(vsite_section,"MOLNAME",i_rep_section=ig,error=error) + END IF + ! Let's tag to which molecule we want to apply constraints + CALL section_vals_val_get(vsite_section,"MOLECULE",i_rep_section=ig,& + n_rep_val=nrep,error=error) + IF (nrep/=0) THEN + CALL section_vals_val_get(vsite_section,"MOLECULE",i_rep_section=ig,& + i_val=cons_info%const_vsite_mol(ig),error=error) + END IF + CALL section_vals_val_get(vsite_section,"MOLNAME",i_rep_section=ig,& + n_rep_val=nrep,error=error) + IF (nrep/=0) THEN + CALL section_vals_val_get(vsite_section,"MOLNAME",i_rep_section=ig,& + c_val=cons_info%const_vsite_molname(ig),error=error) + END IF + IF ((cons_info%const_vsite_mol(ig)/=0).AND.(cons_info%const_vsite_molname(ig)/="UNDEF")) THEN + CPPostcondition(.FALSE.,cp_failure_level,routineP,error,failure) + END IF + IF ((cons_info%const_vsite_mol(ig)==0).AND.(cons_info%const_vsite_molname(ig)=="UNDEF").AND.& + (.NOT.cons_info%vsite_intermolecular(ig))) THEN + CPPostcondition(.FALSE.,cp_failure_level,routineP,error,failure) + END IF + CALL section_vals_val_get(vsite_section,"ATOMS",i_rep_section=ig,& + i_vals=ilist,error=error) + CALL section_vals_val_get(vsite_section,"PARAMETERS",i_rep_section=ig,& + r_vals=rlist,error=error) + cons_info%const_vsite_a(ig) = ilist(1) + cons_info%const_vsite_b(ig) = ilist(2) + cons_info%const_vsite_c(ig) = ilist(3) + cons_info%const_vsite_d(ig) = ilist(4) + cons_info%const_vsite_wbc(ig) = rlist(1) + cons_info%const_vsite_wdc(ig) = rlist(2) + END DO + END IF ! FIXED ATOMS CALL section_vals_get(fix_atom_section,explicit=explicit, n_repetition=ncons, error=error) IF (explicit) THEN diff --git a/src/topology_types.F b/src/topology_types.F index 17976fe887..5246f4f02d 100644 --- a/src/topology_types.F +++ b/src/topology_types.F @@ -9,6 +9,8 @@ !> change name from input_file_name to coord_... !> added topology file !> added atom_names +!> Teodoro Laino [tlaino] 12.2008 - Preparing for VIRTUAL SITE constraints +!> (patch by Marcel Baer) !> \author CJM & JGH ! ***************************************************************************** MODULE topology_types @@ -119,6 +121,20 @@ MODULE topology_types LOGICAL, POINTER :: g46_restraint(:) ! Restraints control REAL(KIND=dp), POINTER :: g46_k0(:) ! Restraints control LOGICAL, POINTER, DIMENSION(:) :: g46_exclude_qm, g46_exclude_mm + ! virtual_site + INTEGER :: nconst_vsite + INTEGER, POINTER :: const_vsite_mol(:) + CHARACTER(LEN=default_string_length), POINTER :: const_vsite_molname(:) + INTEGER, POINTER :: const_vsite_a(:) + INTEGER, POINTER :: const_vsite_b(:) + INTEGER, POINTER :: const_vsite_c(:) + INTEGER, POINTER :: const_vsite_d(:) + REAL(KIND=dp), POINTER :: const_vsite_wbc(:) + REAL(KIND=dp), POINTER :: const_vsite_wdc(:) + LOGICAL, POINTER :: vsite_intermolecular(:) + LOGICAL, POINTER :: vsite_restraint(:) ! Restraints control + REAL(KIND=dp), POINTER :: vsite_k0(:) ! Restraints control + LOGICAL, POINTER, DIMENSION(:) :: vsite_exclude_qm, vsite_exclude_mm END TYPE constraint_info_type ! ***************************************************************************** @@ -151,6 +167,7 @@ MODULE topology_types LOGICAL :: const_colv LOGICAL :: const_33 LOGICAL :: const_46 + LOGICAL :: const_vsite LOGICAL :: charge_occup LOGICAL :: charge_beta LOGICAL :: charge_extended @@ -226,6 +243,7 @@ CONTAINS topology%const_colv=.FALSE. topology%const_33=.FALSE. topology%const_46=.FALSE. + topology%const_vsite=.FALSE. topology%charge_occup=.FALSE. topology%charge_beta=.FALSE. topology%charge_extended=.FALSE. @@ -360,6 +378,20 @@ CONTAINS NULLIFY(constraint_info%g46_restraint) NULLIFY(constraint_info%g46_k0) NULLIFY(constraint_info%g46_exclude_qm, constraint_info%g46_exclude_mm) + ! virtual_site + constraint_info%nconst_vsite = 0 + NULLIFY(constraint_info%const_vsite_mol) + NULLIFY(constraint_info%const_vsite_molname) + NULLIFY(constraint_info%const_vsite_a) + NULLIFY(constraint_info%const_vsite_b) + NULLIFY(constraint_info%const_vsite_c) + NULLIFY(constraint_info%const_vsite_d) + NULLIFY(constraint_info%const_vsite_wbc) + NULLIFY(constraint_info%const_vsite_wdc) + NULLIFY(constraint_info%vsite_intermolecular) + NULLIFY(constraint_info%vsite_restraint) + NULLIFY(constraint_info%vsite_k0) + NULLIFY(constraint_info%vsite_exclude_qm, constraint_info%vsite_exclude_mm) END SUBROUTINE init_constraint @@ -798,6 +830,59 @@ CONTAINS DEALLOCATE(constraint_info%g46_exclude_mm, stat=stat) CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) END IF + ! virtual_site + IF(ASSOCIATED(constraint_info%const_vsite_mol)) THEN + DEALLOCATE(constraint_info%const_vsite_mol,STAT=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + END IF + IF(ASSOCIATED(constraint_info%const_vsite_molname)) THEN + DEALLOCATE(constraint_info%const_vsite_molname,STAT=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + END IF + IF(ASSOCIATED(constraint_info%const_vsite_a)) THEN + DEALLOCATE(constraint_info%const_vsite_a,STAT=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + END IF + IF(ASSOCIATED(constraint_info%const_vsite_b)) THEN + DEALLOCATE(constraint_info%const_vsite_b,STAT=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + END IF + IF(ASSOCIATED(constraint_info%const_vsite_c)) THEN + DEALLOCATE(constraint_info%const_vsite_c,STAT=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + END IF + IF(ASSOCIATED(constraint_info%const_vsite_d)) THEN + DEALLOCATE(constraint_info%const_vsite_d,STAT=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + END IF + IF(ASSOCIATED(constraint_info%const_vsite_wbc)) THEN + DEALLOCATE(constraint_info%const_vsite_wbc,STAT=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + END IF + IF(ASSOCIATED(constraint_info%const_vsite_wdc)) THEN + DEALLOCATE(constraint_info%const_vsite_wdc,STAT=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + END IF + IF(ASSOCIATED(constraint_info%vsite_intermolecular)) THEN + DEALLOCATE(constraint_info%vsite_intermolecular,STAT=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + END IF + IF(ASSOCIATED(constraint_info%vsite_restraint)) THEN + DEALLOCATE(constraint_info%vsite_restraint,STAT=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + END IF + IF(ASSOCIATED(constraint_info%vsite_k0)) THEN + DEALLOCATE(constraint_info%vsite_k0,STAT=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + END IF + IF(ASSOCIATED(constraint_info%vsite_exclude_qm)) THEN + DEALLOCATE(constraint_info%vsite_exclude_qm, stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + END IF + IF(ASSOCIATED(constraint_info%vsite_exclude_mm)) THEN + DEALLOCATE(constraint_info%vsite_exclude_mm, stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + END IF END IF END SUBROUTINE deallocate_constraint diff --git a/tests/Fist/regtest-10/TEST_FILES_RESET b/tests/Fist/regtest-10/TEST_FILES_RESET index cca4f02a82..2a7f3d87f9 100644 --- a/tests/Fist/regtest-10/TEST_FILES_RESET +++ b/tests/Fist/regtest-10/TEST_FILES_RESET @@ -31,3 +31,7 @@ nh3-meta-2.inp # numerical instabilities nh3-meta-3.inp nh3-meta-4.inp +# numerics +nh3-meta-2.inp +nh3-meta-3.inp +nh3-meta-4.inp diff --git a/tests/Fist/regtest-4/TEST_FILES_RESET b/tests/Fist/regtest-4/TEST_FILES_RESET index cce41647c1..0d529fcf93 100644 --- a/tests/Fist/regtest-4/TEST_FILES_RESET +++ b/tests/Fist/regtest-4/TEST_FILES_RESET @@ -177,3 +177,8 @@ nh3-meta-walks_2.inp nh3-meta-walks_1r.inp # unified cube_center nh3-meta-walks_2r.inp +# numerics +nh3-meta-walks_1.inp +nh3-meta-walks_1r.inp +nh3-meta-walks_2.inp +nh3-meta-walks_2r.inp diff --git a/tests/Fist/regtest-5/TEST_FILES b/tests/Fist/regtest-5/TEST_FILES index 7b52a7905e..2b1bf1cdca 100644 --- a/tests/Fist/regtest-5/TEST_FILES +++ b/tests/Fist/regtest-5/TEST_FILES @@ -40,3 +40,5 @@ H2OFe-5_75.inp 2 acn-2.inp 0 H2OFe-5_75-2.inp 0 I3.inp 2 +#TIP4P +water_tip4p.inp 2 diff --git a/tests/Fist/regtest-5/water_tip4p.inp b/tests/Fist/regtest-5/water_tip4p.inp new file mode 100644 index 0000000000..7d30e8fbab --- /dev/null +++ b/tests/Fist/regtest-5/water_tip4p.inp @@ -0,0 +1,224 @@ +&FORCE_EVAL + METHOD FIST + &MM + &FORCEFIELD + &SPLINE + EMAX_SPLINE 1000000.0 + RCUT_NB 12.0 + &END SPLINE +!GROMACS parms +![ atomtypes ] +!; full atom descriptions are available in ffoplsaa.atp +!; name bond_type mass charge ptype sigma epsilon +! opls_113 OW 15.99940 0.000 A 3.15365e-01 6.48520e-01 +! opls_114 HW 1.00800 0.520 A 0.00000e+00 0.00000e+00 +! opls_115 MW 0.00000 -1.040 D 0.00000e+00 0.00000e+00 +! +![ bondtypes ] +!; i j func b0 kb +! OW HW 1 0.09572 502080.0 ; For TIP4F Water - wlj 1/98 +! OW LP 1 0.01750 753120.0 ; -idem- +! +![ angletypes ] +!; i j k func th0 cth +! HW OW HW 1 109.500 627.600 ; For TIP4F Water - wj 1/98 +! HW OW LP 1 54.750 418.400 ; For TIP4F Water - wj 1/98 + &BEND + ATOMS HW OW HW + THETA0 [deg] 109.500 + !K [rad^2kjmol] 627.600 + K [rad^2kjmol] 0.000 + KIND G87 + &END BEND + &BEND + ATOMS HW OW MW + THETA0 [deg] 54.750 + !K [rad^2kjmol] 418.400 + K [rad^2kjmol] 0.000 + KIND G87 + &END BEND + &BEND + ATOMS MW HW OW + THETA0 [deg] 109.500 + K [rad^2kjmol] 0.000 + KIND G87 + &END BEND + &BEND + ATOMS MW HW HW + THETA0 [deg] 109.500 + K [rad^2kjmol] 0.000 + KIND G87 + &END BEND + &BEND + ATOMS OW HW HW + THETA0 [deg] 109.500 + K [rad^2kjmol] 0.000 + KIND G87 + &END BEND + &BEND + ATOMS OW MW HW + THETA0 [deg] 109.500 + K [rad^2kjmol] 0.000 + KIND G87 + &END BEND + &BEND + ATOMS HW MW HW + THETA0 [deg] 109.500 + K [rad^2kjmol] 0.000 + KIND G87 + &END BEND + &BOND + ATOMS OW HW + K [nm^-2kjmol] 502080.0 + !K [nm^-2kjmol] 0.0 + R0 [nm] 0.09572 + KIND G87 + &END BOND + &BOND + ATOMS OW MW + !K [nm^-2kjmol] 753120.0 + K [nm^-2kjmol] 00.0 + R0 [nm] 0.01750 + KIND G87 + &END BOND + &BOND + ATOMS HW HW + K [nm^-2kjmol] 000000.0 + R0 [nm] 0.01750 + KIND G87 + &END BOND + &BOND + ATOMS HW MW + K [nm^-2kjmol] 000000.0 + R0 [nm] 0.01750 + KIND G87 + &END BOND + &NONBONDED + &LENNARD-JONES + ATOMS OW OW + EPSILON [kjmol] 6.36386E-01 + SIGMA [nm] 3.15061E-01 + &END LENNARD-JONES + &LENNARD-JONES + ATOMS OW MW + EPSILON 0.0 + SIGMA 0.1 + &END LENNARD-JONES + &LENNARD-JONES + ATOMS OW HW + EPSILON 0.0 + SIGMA 0.1 + &END LENNARD-JONES + &LENNARD-JONES + ATOMS HW HW + EPSILON 0.0 + SIGMA 0.1 + &END LENNARD-JONES + &LENNARD-JONES + ATOMS HW MW + EPSILON 0.0 + SIGMA 0.1 + &END LENNARD-JONES + &LENNARD-JONES + ATOMS MW MW + EPSILON 0.0 + SIGMA 0.1 + &END LENNARD-JONES + &END NONBONDED + &CHARGE + ATOM OW + ! CHARGE -1.04 + CHARGE 0.000000 + &END CHARGE + &CHARGE + ATOM HW + CHARGE 0.52 + !CHARGE 0.00 + &END CHARGE + &CHARGE + ATOM MW + CHARGE -1.04 + ! CHARGE 0.00 + &END CHARGE + &END FORCEFIELD + &POISSON + &EWALD + !EWALD_TYPE spme + EWALD_TYPE NONE + ALPHA .44 + GMAX 24 + O_SPLINE 6 + &END EWALD + &END POISSON + &END MM + &SUBSYS + &CELL + ABC 40.000 40.000 40.000 + &END CELL + &TOPOLOGY + COORD_FILE_NAME ../sample_pdb/water_tip4p.pdb + COORDINATE pdb + CONN_FILE_NAME ../sample_psf/water_tip4p.psf + CONN_FILE_FORMAT PSF + &END TOPOLOGY + &KIND OW + ELEMENT O + MASS 15.99940 + &END KIND + &KIND HW + ELEMENT H + MASS 1.00800 + &END KIND + &KIND MW + ELEMENT H + MASS 0.00000 + &END KIND + &PRINT + &ATOMIC_COORDINATES SILENT + &END ATOMIC_COORDINATES + &MOLECULES ON + &END MOLECULES + &END PRINT + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT tip4p + RUN_TYPE MD +&END GLOBAL +&MOTION + &MD + ENSEMBLE NVE + STEPS 200 + TIMESTEP 1.0 + TEMPERATURE 100 + &END MD + &CONSTRAINT + &G3X3 + INTERMOLECULAR FALSE + DISTANCES 1.89035 1.89035 3.0869 + ATOMS 1 2 3 + MOLECULE 1 + &END G3X3 + &VIRTUAL_SITE + INTERMOLECULAR FALSE + ATOMS 4 2 1 3 + PARAMETERS 0.128012065 0.128012065 + MOLECULE 1 + &END VIRTUAL_SITE + &END CONSTRAINT + &PRINT + &FORCES ON + &EACH + MD 2 + &END EACH + FORMAT XYZ + &END FORCES + &VELOCITIES ON + &EACH + MD 2 + &END EACH + FORMAT XYZ + &END VELOCITIES + &END PRINT +&END MOTION + diff --git a/tests/Fist/sample_pdb/water_tip4p.pdb b/tests/Fist/sample_pdb/water_tip4p.pdb new file mode 100644 index 0000000000..06d8f478c5 --- /dev/null +++ b/tests/Fist/sample_pdb/water_tip4p.pdb @@ -0,0 +1,10 @@ +CRYST1 0.000 0.000 0.000 90.00 90.00 90.00 P 1 1 +ATOM 1 OW TIP4X 0 4.257 3.685 7.048 1.00 0.00 O +ATOM 2 HW TIP4X 0 3.505 3.065 6.819 1.00 0.00 H +ATOM 3 HW TIP4X 0 4.164 4.529 6.518 1.00 0.00 H +ATOM 4 MW TIP4X 0 4.149 3.714 6.951 1.00 0.00 H +ATOM 5 OW TIP4X 1 6.650 2.043 6.264 1.00 0.00 O +ATOM 6 HW TIP4X 1 5.872 2.583 6.586 1.00 0.00 H +ATOM 7 HW TIP4X 1 7.297 1.910 7.015 1.00 0.00 H +ATOM 8 MW TIP4X 1 6.633 2.095 6.401 1.00 0.00 H +END diff --git a/tests/Fist/sample_psf/water_tip4p.psf b/tests/Fist/sample_psf/water_tip4p.psf new file mode 100644 index 0000000000..7e503592b5 --- /dev/null +++ b/tests/Fist/sample_psf/water_tip4p.psf @@ -0,0 +1,41 @@ +PSF EXT + + 1 !NTITLE + CP2K generated DUMP of connectivity + + 8 !NATOM + 1 MOL1 1 TIP OW OW 0.000000 15.999 0 + 2 MOL1 1 TIP HW HW 0.520000 1.008 0 + 3 MOL1 1 TIP HW HW 0.520000 1.008 0 + 4 MOL1 1 TIP MW MW -1.040000 0.000 0 + 5 MOL1 2 TIP OW OW 0.000000 15.999 0 + 6 MOL1 2 TIP HW HW 0.520000 1.008 0 + 7 MOL1 2 TIP HW HW 0.520000 1.008 0 + 8 MOL1 2 TIP MW MW -1.040000 0.000 0 + + + 10 !NBOND + 1 2 1 3 5 6 5 7 + 4 2 4 3 8 6 8 7 + 4 1 8 5 + + + 2 !NTHETA + 2 1 3 6 5 7 + + + 0 !NPHI + + + 0 !NIMPHI + + + 0 !NDON + + + 0 !NACC + + + 0 !NNB + + diff --git a/tests/SE/regtest-2/TEST_FILES_RESET b/tests/SE/regtest-2/TEST_FILES_RESET index bac0302897..39bacfc4cc 100644 --- a/tests/SE/regtest-2/TEST_FILES_RESET +++ b/tests/SE/regtest-2/TEST_FILES_RESET @@ -139,3 +139,5 @@ md.inp scan_zn_1w_cp2k.inp # numerical instabilities hcl_ot.inp +# numerics +hcl_ot.inp