diff --git a/src/input_cp2k_motion.F b/src/input_cp2k_motion.F index c734a6b909..06a450400f 100644 --- a/src/input_cp2k_motion.F +++ b/src/input_cp2k_motion.F @@ -900,6 +900,20 @@ CONTAINS CALL section_add_keyword(section,keyword,error=error) CALL keyword_release(keyword,error=error) + 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_collective_section @@ -993,6 +1007,22 @@ CONTAINS CALL section_add_keyword(section,keyword,error=error) CALL keyword_release(keyword,error=error) + CALL keyword_create(keyword, name="EXCLUDE_QM",& + description="Does not apply the constraint to the QM region within a QM/MM calculation."//& + " This keyword is active only together with MOLNAME",& + 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."//& + " This keyword is active only together with MOLNAME",& + 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_fixed_atom_section @@ -1070,10 +1100,25 @@ CONTAINS ! Real CALL keyword_create(keyword, name="DISTANCES",& description="The constrained distances' values.",& - usage="DISTANCES {real} {real} {real}", type_of_var=real_t, n_var=-1,error=error) + usage="DISTANCES {real} {real} {real}", type_of_var=real_t,& + unit_str="internal_cp2k",n_var=-1,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_g3x3_section @@ -1244,10 +1289,24 @@ CONTAINS CALL keyword_create(keyword, name="DISTANCES",& description="The constrained distances' values.",& usage="DISTANCES {real} {real} {real} {real} {real} {real}",& - type_of_var=real_t, n_var=6,error=error) + type_of_var=real_t, n_var=6,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_g4x6_section diff --git a/src/md.F b/src/md.F index 933068a0b2..5c789e269a 100644 --- a/src/md.F +++ b/src/md.F @@ -159,6 +159,9 @@ CONTAINS TYPE(section_vals_type), POINTER :: root_section TYPE(cp_error_type), INTENT(inout) :: error + CHARACTER(len=*), PARAMETER :: routineN = 'read_md_section', & + routineP = moduleN//':'//routineN + CHARACTER(LEN=default_path_length) :: filename INTEGER :: iw INTEGER, DIMENSION(:), POINTER :: iprint @@ -313,6 +316,14 @@ CONTAINS ' SHOCK| Barostat Mass [', 'au', '] ', & simpar % cmass END IF + ! Print warning for temp_tol + IF (simpar % temp_tol > 0.0_dp) THEN + CALL cp_assert(.FALSE.,cp_warning_level,cp_assertion_failed,routineP,& + "A temperature tolerance (TEMP_TOL) is used during the MD. "//& + "Due to the velocity rescaling algorithm jumps may appear in the conserved quantity. "//& +CPSourceFileRef,& + only_ionode=.TRUE.) + END IF END IF CALL cp_print_key_finished_output(iw,logger,md_section,& diff --git a/src/topology_constraint_util.F b/src/topology_constraint_util.F index 10cce5d0ed..76390dd23a 100644 --- a/src/topology_constraint_util.F +++ b/src/topology_constraint_util.F @@ -132,8 +132,9 @@ CONTAINS INTEGER, DIMENSION(:), POINTER :: constr_x_glob, inds, & molecule_list LOGICAL :: exclude_mm, exclude_qm, failure, fix_atom_mm, & - fix_atom_molname, fix_atom_qm, fix_atom_qmmm, fix_fixed_atom, ionode, & - is_qm, ishbond, ldummy, restart_restraint_clv, restart_restraint_pos + fix_atom_molname, fix_atom_qm, fix_atom_qmmm, fix_fixed_atom, & + found_molname, ionode, is_qm, ishbond, ldummy, restart_restraint_clv, & + restart_restraint_pos LOGICAL, ALLOCATABLE, DIMENSION(:) :: missed_molname REAL(KIND=dp) :: rmod, rvec(3) REAL(KIND=dp), DIMENSION(:), POINTER :: hdist, r @@ -252,12 +253,17 @@ CONTAINS END DO CALL section_vals_val_get(hbonds_section,"MOLNAME",c_vals=cnds,error=error) DO i = 1, SIZE(cnds) + found_molname = .FALSE. DO k = 1, SIZE(molecule_kind_set) molecule_kind => molecule_kind_set(k) name = molecule_kind%name ldummy = qmmm_ff_precond_only_qm(id1=name) - IF (cnds(i) == name) constr_x_mol(k)%constr(1) = 1 + IF (cnds(i) == name) THEN + constr_x_mol(k)%constr(1) = 1 + found_molname = .TRUE. + END IF END DO + CALL print_warning_molname(found_molname, cnds(i), error) END DO END IF END IF @@ -321,8 +327,12 @@ CONTAINS CALL reallocate(cons_info%colv_intermolecular,1,cons_info%nconst_colv) CALL reallocate(cons_info%colv_restraint,1,cons_info%nconst_colv) CALL reallocate(cons_info%colv_k0,1,cons_info%nconst_colv) + CALL reallocate(cons_info%colv_exclude_qm,1,cons_info%nconst_colv) + CALL reallocate(cons_info%colv_exclude_mm,1,cons_info%nconst_colv) ! Bonds involving hydrogens are by their nature only intramolecular cons_info%colv_intermolecular(n_start_colv+1:cons_info%nconst_colv) = .FALSE. + cons_info%colv_exclude_qm(n_start_colv+1:cons_info%nconst_colv) = .FALSE. + cons_info%colv_exclude_mm(n_start_colv+1:cons_info%nconst_colv) = .FALSE. cons_info%colv_restraint(n_start_colv+1:cons_info%nconst_colv) = cons_info%hbonds_restraint cons_info%colv_k0(n_start_colv+1:cons_info%nconst_colv) = cons_info%hbonds_k0 ! @@ -414,6 +424,8 @@ CONTAINS constr_x_mol,& constr_x_glob,& molecule_kind_set,& + cons_info%colv_exclude_qm,& + cons_info%colv_exclude_mm,& error) ! Intramolecular constraints gind = 0 @@ -485,6 +497,8 @@ CONTAINS constr_x_mol,& constr_x_glob,& molecule_kind_set,& + cons_info%g33_exclude_qm,& + cons_info%g33_exclude_mm,& error) ! Intramolecular constraints DO ii = 1, SIZE(molecule_kind_set) @@ -542,6 +556,8 @@ CONTAINS constr_x_mol,& constr_x_glob,& molecule_kind_set,& + cons_info%g46_exclude_qm,& + cons_info%g46_exclude_mm,& error) ! Intramolecular constraints DO ii = 1, SIZE(molecule_kind_set) @@ -604,7 +620,7 @@ CONTAINS molecule_kind => molecule_kind_set(i) CALL get_molecule_kind(molecule_kind=molecule_kind,& nmolecule=nmolecule,molecule_list=molecule_list, name=molname) - ldummy = qmmm_ff_precond_only_qm(id1=molname) + is_qm = qmmm_ff_precond_only_qm(id1=molname) WHERE(molname.EQ.cons_info%fixed_molnames) missed_molname=.FALSE. END WHERE @@ -615,7 +631,13 @@ CONTAINS CALL get_molecule ( molecule, first_atom = first, last_atom = last ) fix_atom_molname=.FALSE. IF (ASSOCIATED(cons_info%fixed_molnames)) THEN - IF (ANY(cons_info%fixed_molnames.EQ.molname)) fix_atom_molname=.TRUE. + DO k = 1, SIZE(cons_info%fixed_molnames) + IF (cons_info%fixed_molnames(k).EQ.molname) THEN + fix_atom_molname=.TRUE. + IF (is_qm.AND.cons_info%fixed_exclude_qm(k)) fix_atom_molname=.FALSE. + IF ((.NOT.is_qm).AND.cons_info%fixed_exclude_mm(k)) fix_atom_molname=.FALSE. + END IF + END DO ENDIF DO k = first, last fix_atom_qmmm=.FALSE. @@ -642,111 +664,113 @@ CONTAINS CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) kk = 0 nfixd_restraint = 0 - DO j=1,nmolecule - molecule => molecule_set(molecule_list(j)) - CALL get_molecule ( molecule, first_atom = first, last_atom = last ) - fix_atom_molname=.FALSE. - IF (ASSOCIATED(cons_info%fixed_molnames)) THEN - DO k1loc = 1, SIZE(cons_info%fixed_molnames) - IF (cons_info%fixed_molnames(k1loc).EQ.molname) THEN - fix_atom_molname=.TRUE. - EXIT + IF (nfixed_atoms/=0) THEN + DO j=1,nmolecule + molecule => molecule_set(molecule_list(j)) + CALL get_molecule ( molecule, first_atom = first, last_atom = last ) + fix_atom_molname=.FALSE. + IF (ASSOCIATED(cons_info%fixed_molnames)) THEN + DO k1loc = 1, SIZE(cons_info%fixed_molnames) + IF (cons_info%fixed_molnames(k1loc).EQ.molname) THEN + fix_atom_molname=.TRUE. + EXIT + END IF + END DO + ENDIF + DO k=first,last + ! FIXED LIST ATOMS + fix_fixed_atom = .FALSE. + DO k2loc=1,SIZE(cons_info%fixed_atoms) + IF (cons_info%fixed_atoms(k2loc) == k) THEN + fix_fixed_atom = .TRUE. + EXIT + END IF + END DO + ! QMMM FIXED ATOMS (QM OR MM) + fix_atom_qmmm=.FALSE. + fix_atom_mm =.FALSE. + fix_atom_qm =.FALSE. + IF (PRESENT(qmmm_env)) THEN + SELECT CASE(cons_info%freeze_qm) + CASE(do_constr_atomic) + IF (ANY(qmmm_env%qm_atom_index==k)) THEN + fix_atom_qmmm=.TRUE. + fix_atom_qm = .TRUE. + END IF + CASE(do_constr_molec) + IF (ANY(qmmm_env%qm_molecule_index == molecule_list(j))) THEN + fix_atom_qmmm = .TRUE. + fix_atom_qm = .TRUE. + END IF + END SELECT + SELECT CASE(cons_info%freeze_mm) + CASE(do_constr_atomic) + IF (ALL(qmmm_env%qm_atom_index/=k)) THEN + fix_atom_qmmm=.TRUE. + fix_atom_mm = .TRUE. + END IF + CASE(do_constr_molec) + IF (ALL(qmmm_env%qm_molecule_index /= molecule_list(j))) THEN + fix_atom_qmmm = .TRUE. + fix_atom_mm = .TRUE. + END IF + END SELECT + ! We should never reach this point but let's check it anyway + IF (fix_atom_qm.AND.fix_atom_mm) THEN + CALL stop_program(TRIM(routineN),"Atom number: "//cp_to_string(k)//& + " has been defined both QM and MM. General Error!") + END IF END IF - END DO - ENDIF - DO k=first,last - ! FIXED LIST ATOMS - fix_fixed_atom = .FALSE. - DO k2loc=1,SIZE(cons_info%fixed_atoms) - IF (cons_info%fixed_atoms(k2loc) == k) THEN - fix_fixed_atom = .TRUE. - EXIT - END IF - END DO - ! QMMM FIXED ATOMS (QM OR MM) - fix_atom_qmmm=.FALSE. - fix_atom_mm =.FALSE. - fix_atom_qm =.FALSE. - IF (PRESENT(qmmm_env)) THEN - SELECT CASE(cons_info%freeze_qm) - CASE(do_constr_atomic) - IF (ANY(qmmm_env%qm_atom_index==k)) THEN - fix_atom_qmmm=.TRUE. - fix_atom_qm = .TRUE. - END IF - CASE(do_constr_molec) - IF (ANY(qmmm_env%qm_molecule_index == molecule_list(j))) THEN - fix_atom_qmmm = .TRUE. - fix_atom_qm = .TRUE. - END IF - END SELECT - SELECT CASE(cons_info%freeze_mm) - CASE(do_constr_atomic) - IF (ALL(qmmm_env%qm_atom_index/=k)) THEN - fix_atom_qmmm=.TRUE. - fix_atom_mm = .TRUE. - END IF - CASE(do_constr_molec) - IF (ALL(qmmm_env%qm_molecule_index /= molecule_list(j))) THEN - fix_atom_qmmm = .TRUE. - fix_atom_mm = .TRUE. - END IF - END SELECT - ! We should never reach this point but let's check it anyway - IF (fix_atom_qm.AND.fix_atom_mm) THEN + ! Check that the fixed atom constraint/restraint is unique + IF ((fix_fixed_atom.AND.fix_atom_qmmm).OR.(fix_fixed_atom.AND.fix_atom_molname)& + .OR.(fix_atom_qmmm.AND.fix_atom_molname)) THEN CALL stop_program(TRIM(routineN),"Atom number: "//cp_to_string(k)//& - " has been defined both QM and MM. General Error!") + " has been constrained/restrained to be fixed in more than one"//& + " input section. Check and correct your input file!") END IF - END IF - ! Check that the fixed atom constraint/restraint is unique - IF ((fix_fixed_atom.AND.fix_atom_qmmm).OR.(fix_fixed_atom.AND.fix_atom_molname)& - .OR.(fix_atom_qmmm.AND.fix_atom_molname)) THEN - CALL stop_program(TRIM(routineN),"Atom number: "//cp_to_string(k)//& - " has been constrained/restrained to be fixed in more than one"//& - " input section. Check and correct your input file!") - END IF - ! Let's store the atom index - IF (fix_fixed_atom.OR.fix_atom_qmmm.OR.fix_atom_molname) THEN - kk = kk + 1 - fixd_list(kk)%fixd = k - fixd_list(kk)%coord = particle_set(k)%r - ! Possibly Restraint - IF (fix_fixed_atom) THEN - fixd_list(kk)%restraint%active = cons_info%fixed_restraint(k2loc) - fixd_list(kk)%restraint%k0 = cons_info%fixed_k0(k2loc) - ELSEIF (fix_atom_qm) THEN - fixd_list(kk)%restraint%active = cons_info%fixed_qm_restraint - fixd_list(kk)%restraint%k0 = cons_info%fixed_qm_k0 - ELSEIF (fix_atom_mm) THEN - fixd_list(kk)%restraint%active = cons_info%fixed_mm_restraint - fixd_list(kk)%restraint%k0 = cons_info%fixed_mm_k0 - ELSEIF (fix_atom_molname) THEN - fixd_list(kk)%restraint%active = cons_info%fixed_mol_restraint(k1loc) - fixd_list(kk)%restraint%k0 = cons_info%fixed_mol_k0(k1loc) - ELSE - ! Should never reach this point - CPPostcondition(.FALSE.,cp_failure_level,routineP,error,failure) - END IF - IF (fixd_list(kk)%restraint%active) THEN - nfixd_restraint = nfixd_restraint+ 1 - nfixd_restart = nfixd_restart + 1 - IF (restart_restraint_pos) THEN - ! Read coord0 value for restraint - CALL section_vals_val_get(fixd_restr_rest,"_DEFAULT_KEYWORD_",& - i_rep_val=nfixd_restart, r_vals=r, error=error) - fixd_list(kk)%coord = r + ! Let's store the atom index + IF (fix_fixed_atom.OR.fix_atom_qmmm.OR.fix_atom_molname) THEN + kk = kk + 1 + fixd_list(kk)%fixd = k + fixd_list(kk)%coord = particle_set(k)%r + ! Possibly Restraint + IF (fix_fixed_atom) THEN + fixd_list(kk)%restraint%active = cons_info%fixed_restraint(k2loc) + fixd_list(kk)%restraint%k0 = cons_info%fixed_k0(k2loc) + ELSEIF (fix_atom_qm) THEN + fixd_list(kk)%restraint%active = cons_info%fixed_qm_restraint + fixd_list(kk)%restraint%k0 = cons_info%fixed_qm_k0 + ELSEIF (fix_atom_mm) THEN + fixd_list(kk)%restraint%active = cons_info%fixed_mm_restraint + fixd_list(kk)%restraint%k0 = cons_info%fixed_mm_k0 + ELSEIF (fix_atom_molname) THEN + fixd_list(kk)%restraint%active = cons_info%fixed_mol_restraint(k1loc) + fixd_list(kk)%restraint%k0 = cons_info%fixed_mol_k0(k1loc) ELSE - ! Write coord0 value for restraint - ALLOCATE(r(3),stat=stat) - CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) - r = fixd_list(kk)%coord - CALL section_vals_val_set(fixd_restr_rest,"_DEFAULT_KEYWORD_",& - i_rep_val=nfixd_restart, r_vals_ptr=r, error=error) + ! Should never reach this point + CPPostcondition(.FALSE.,cp_failure_level,routineP,error,failure) + END IF + IF (fixd_list(kk)%restraint%active) THEN + nfixd_restraint = nfixd_restraint+ 1 + nfixd_restart = nfixd_restart + 1 + IF (restart_restraint_pos) THEN + ! Read coord0 value for restraint + CALL section_vals_val_get(fixd_restr_rest,"_DEFAULT_KEYWORD_",& + i_rep_val=nfixd_restart, r_vals=r, error=error) + fixd_list(kk)%coord = r + ELSE + ! Write coord0 value for restraint + ALLOCATE(r(3),stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + r = fixd_list(kk)%coord + CALL section_vals_val_set(fixd_restr_rest,"_DEFAULT_KEYWORD_",& + i_rep_val=nfixd_restart, r_vals_ptr=r, error=error) + END IF END IF END IF - END IF + END DO END DO - END DO + END IF IF (iw>0) THEN WRITE(iw,*)"MOLECULE KIND:",i," NR. FIXED ATOMS:",SIZE(fixd_list(:)%fixd)," LIST::",fixd_list(:)%fixd END IF @@ -1225,7 +1249,7 @@ END SUBROUTINE setup_lg4x6 !! INPUTS !!*** ************************************************************************* SUBROUTINE give_constraint_array(const_mol, const_molname, const_intermolecular,& - constr_x_mol, constr_x_glob, molecule_kind_set, error) + constr_x_mol, constr_x_glob, molecule_kind_set, exclude_qm, exclude_mm, error) INTEGER, DIMENSION(:), POINTER :: const_mol CHARACTER(LEN=default_string_length), & @@ -1236,6 +1260,7 @@ END SUBROUTINE setup_lg4x6 INTEGER, DIMENSION(:), POINTER :: constr_x_glob TYPE(molecule_kind_type), DIMENSION(:), & POINTER :: molecule_kind_set + LOGICAL, DIMENSION(:), POINTER :: exclude_qm, exclude_mm TYPE(cp_error_type), INTENT(inout) :: error CHARACTER(len=*), PARAMETER :: routineN = 'give_constraint_array', & @@ -1244,7 +1269,7 @@ END SUBROUTINE setup_lg4x6 CHARACTER(LEN=default_string_length) :: myname, name INTEGER :: handle, i, iglob, isize, k, & stat - LOGICAL :: failure, ldummy + LOGICAL :: failure, found_molname, is_qm TYPE(molecule_kind_type), POINTER :: molecule_kind failure = .FALSE. @@ -1275,16 +1300,21 @@ END SUBROUTINE setup_lg4x6 constr_x_mol(k)%constr(isize+1) = i ELSE myname = const_molname(i) + found_molname = .FALSE. DO k = 1, SIZE(molecule_kind_set) molecule_kind => molecule_kind_set(k) name = molecule_kind%name - ldummy = qmmm_ff_precond_only_qm(id1=name) + is_qm = qmmm_ff_precond_only_qm(id1=name) + IF (is_qm.AND.exclude_qm(i)) CYCLE + IF (.NOT.is_qm.AND.exclude_mm(i)) CYCLE IF (name==myname) THEN isize = SIZE(constr_x_mol(k)%constr) CALL reallocate(constr_x_mol(k)%constr, 1, isize+1) constr_x_mol(k)%constr(isize+1) = i + found_molname = .TRUE. END IF END DO + CALL print_warning_molname(found_molname, myname, error) END IF END IF END DO @@ -1292,4 +1322,47 @@ END SUBROUTINE setup_lg4x6 CALL timestop(0.0_dp,handle) END SUBROUTINE give_constraint_array +!****************************************************************************** +!!****** topology_constraint_util/print_warning_molname [1.0] * +!! +!! NAME +!! print_warning_molname +!! +!! SYNOPSIS +!! +!! FUNCTION +!! Prints a warning message if undefined molnames are used to define constraints +!! +!! AUTHOR +!! Teodoro Laino [2007] - Zurich University +!! +!! MODIFICATION HISTORY +!! +!! +!! OPTIONS +!! +!! INPUTS +!!*** ************************************************************************* + SUBROUTINE print_warning_molname(found, name, error) + LOGICAL, INTENT(IN) :: found + CHARACTER(LEN=*), INTENT(IN) :: name + TYPE(cp_error_type), INTENT(inout) :: error + + CHARACTER(len=*), PARAMETER :: routineN = 'print_warning_molname', & + routineP = moduleN//':'//routineN + + LOGICAL :: failure + + failure = .FALSE. + CALL cp_assert(found,cp_warning_level,cp_assertion_failed,routineP,& + " MOLNAME ("//TRIM(name)//") was defined for constraints, but this molecule name "//& + "is not defined. Please check carefully your PDB, PSF (has priority over PDB) or "//& + "input driven CP2K coordinates. In case you may not find the reason for this warning "//& + "it may be a good idea to print all molecule information (including kind name) activating "//& + "the print_key MOLECULES specific of the SUBSYS%PRINT section. "//& +CPSourceFileRef,& + only_ionode=.TRUE.) + + END SUBROUTINE print_warning_molname + END MODULE topology_constraint_util diff --git a/src/topology_input.F b/src/topology_input.F index dcc268c730..8680663e38 100644 --- a/src/topology_input.F +++ b/src/topology_input.F @@ -230,6 +230,10 @@ SUBROUTINE read_constraints_section(topology, colvar_p, constraint_section, erro CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) ALLOCATE(cons_info%g33_k0(ncons),stat=stat) CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + ALLOCATE(cons_info%g33_exclude_qm(ncons),stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + ALLOCATE(cons_info%g33_exclude_mm(ncons),stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) DO ig = 1, ncons CALL check_restraint(g3x3_section,& is_restraint=cons_info%g33_restraint(ig), & @@ -239,6 +243,11 @@ SUBROUTINE read_constraints_section(topology, colvar_p, constraint_section, erro error=error) cons_info%const_g33_mol(ig) = 0 cons_info%const_g33_molname(ig) = "UNDEF" + ! Exclude QM or MM + CALL section_vals_val_get(g3x3_section,"EXCLUDE_QM",i_rep_section=ig,& + l_val=cons_info%g33_exclude_qm(ig),error=error) + CALL section_vals_val_get(g3x3_section,"EXCLUDE_MM",i_rep_section=ig,& + l_val=cons_info%g33_exclude_mm(ig),error=error) ! Intramolecular restraint CALL section_vals_val_get(g3x3_section,"INTERMOLECULAR",i_rep_section=ig,& l_val=cons_info%g33_intermolecular(ig),error=error) @@ -317,6 +326,10 @@ SUBROUTINE read_constraints_section(topology, colvar_p, constraint_section, erro CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) ALLOCATE(cons_info%g46_k0(ncons),stat=stat) CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + ALLOCATE(cons_info%g46_exclude_qm(ncons),stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + ALLOCATE(cons_info%g46_exclude_mm(ncons),stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) DO ig = 1, ncons CALL check_restraint(g4x6_section,& is_restraint=cons_info%g46_restraint(ig), & @@ -326,6 +339,11 @@ SUBROUTINE read_constraints_section(topology, colvar_p, constraint_section, erro error=error) cons_info%const_g46_mol(ig) = 0 cons_info%const_g46_molname(ig) = "UNDEF" + ! Exclude QM or MM + CALL section_vals_val_get(g4x6_section,"EXCLUDE_QM",i_rep_section=ig,& + l_val=cons_info%g46_exclude_qm(ig),error=error) + CALL section_vals_val_get(g4x6_section,"EXCLUDE_MM",i_rep_section=ig,& + l_val=cons_info%g46_exclude_mm(ig),error=error) ! Intramolecular restraint CALL section_vals_val_get(g4x6_section,"INTERMOLECULAR",i_rep_section=ig,& l_val=cons_info%g46_intermolecular(ig),error=error) @@ -389,6 +407,10 @@ SUBROUTINE read_constraints_section(topology, colvar_p, constraint_section, erro CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) ALLOCATE(cons_info%fixed_mol_k0(msize),stat=stat) CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + ALLOCATE(cons_info%fixed_exclude_qm(ncons),stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + ALLOCATE(cons_info%fixed_exclude_mm(ncons),stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) DO ig = 1, ncons isize_old = isize msize_old = msize @@ -439,6 +461,15 @@ SUBROUTINE read_constraints_section(topology, colvar_p, constraint_section, erro CALL reallocate(cons_info%fixed_mol_k0,1,msize+SIZE(tmpstringlist,1)) msize = SIZE(cons_info%fixed_molnames) cons_info%fixed_molnames(:)=tmpstringlist(:) + ! Exclude QM or MM work only if defined MOLNAME + CALL reallocate(cons_info%fixed_exclude_qm,1,msize) + CALL reallocate(cons_info%fixed_exclude_mm,1,msize) + CALL section_vals_val_get(fix_atom_section,"EXCLUDE_QM",i_rep_section=ig,& + l_val=cons_info%fixed_exclude_qm(msize_old+1),error=error) + CALL section_vals_val_get(fix_atom_section,"EXCLUDE_MM",i_rep_section=ig,& + l_val=cons_info%fixed_exclude_mm(msize_old+1),error=error) + cons_info%fixed_exclude_qm(msize_old+1:msize) = cons_info%fixed_exclude_qm(msize_old+1) + cons_info%fixed_exclude_mm(msize_old+1:msize) = cons_info%fixed_exclude_mm(msize_old+1) END IF !Check for restraints IF (n_rep /=0) THEN @@ -513,6 +544,10 @@ SUBROUTINE read_constraints_section(topology, colvar_p, constraint_section, erro CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) ALLOCATE(cons_info%colv_k0(ncons),stat=stat) CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + ALLOCATE(cons_info%colv_exclude_qm(ncons),stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + ALLOCATE(cons_info%colv_exclude_mm(ncons),stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) DO ig = 1, ncons CALL check_restraint(collective_section,& is_restraint=cons_info%colv_restraint(ig), & @@ -522,6 +557,11 @@ SUBROUTINE read_constraints_section(topology, colvar_p, constraint_section, erro error=error) cons_info%const_colv_mol(ig) = 0 cons_info%const_colv_molname(ig) = "UNDEF" + ! Exclude QM or MM + CALL section_vals_val_get(collective_section,"EXCLUDE_QM",i_rep_section=ig,& + l_val=cons_info%colv_exclude_qm(ig),error=error) + CALL section_vals_val_get(collective_section,"EXCLUDE_MM",i_rep_section=ig,& + l_val=cons_info%colv_exclude_mm(ig),error=error) ! Intramolecular restraint CALL section_vals_val_get(collective_section,"INTERMOLECULAR",i_rep_section=ig,& l_val=cons_info%colv_intermolecular(ig),error=error) diff --git a/src/topology_types.F b/src/topology_types.F index 0d4652a4b4..e1d67e45ce 100644 --- a/src/topology_types.F +++ b/src/topology_types.F @@ -85,6 +85,7 @@ MODULE topology_types LOGICAL, POINTER :: fixed_mol_restraint(:) ! Restraints control REAL(KIND=dp), POINTER :: fixed_mol_k0(:) ! Restraints control CHARACTER(LEN=default_string_length), POINTER :: fixed_molnames(:) + LOGICAL, POINTER, DIMENSION(:) :: fixed_exclude_qm, fixed_exclude_mm ! Collective constraints INTEGER :: nconst_colv INTEGER, POINTER :: const_colv_mol(:) @@ -94,6 +95,7 @@ MODULE topology_types LOGICAL, POINTER :: colv_intermolecular(:) LOGICAL, POINTER :: colv_restraint(:) ! Restraints control REAL(KIND=dp), POINTER :: colv_k0(:) ! Restraints control + LOGICAL, POINTER, DIMENSION(:) :: colv_exclude_qm, colv_exclude_mm ! G3x3 INTEGER :: nconst_g33 INTEGER, POINTER :: const_g33_mol(:) @@ -107,6 +109,7 @@ MODULE topology_types LOGICAL, POINTER :: g33_intermolecular(:) LOGICAL, POINTER :: g33_restraint(:) ! Restraints control REAL(KIND=dp), POINTER :: g33_k0(:) ! Restraints control + LOGICAL, POINTER, DIMENSION(:) :: g33_exclude_qm, g33_exclude_mm ! G4x6 INTEGER :: nconst_g46 INTEGER, POINTER :: const_g46_mol(:) @@ -124,6 +127,7 @@ MODULE topology_types LOGICAL, POINTER :: g46_intermolecular(:) LOGICAL, POINTER :: g46_restraint(:) ! Restraints control REAL(KIND=dp), POINTER :: g46_k0(:) ! Restraints control + LOGICAL, POINTER, DIMENSION(:) :: g46_exclude_qm, g46_exclude_mm END TYPE constraint_info_type TYPE topology_parameters_type @@ -165,13 +169,13 @@ MODULE topology_types END TYPE constr_list_type PUBLIC :: atom_info_type,& - connectivity_info_type,& - constraint_info_type,& - topology_parameters_type,& - constr_list_type + connectivity_info_type,& + constraint_info_type,& + topology_parameters_type,& + constr_list_type PUBLIC :: init_topology,& - deallocate_topology + deallocate_topology CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_types' PRIVATE @@ -312,6 +316,7 @@ CONTAINS NULLIFY(constraint_info%fixed_k0) NULLIFY(constraint_info%fixed_mol_restraint) NULLIFY(constraint_info%fixed_mol_k0) + NULLIFY(constraint_info%fixed_exclude_qm, constraint_info%fixed_exclude_mm) ! Collective Constraints constraint_info%nconst_colv = 0 NULLIFY(constraint_info%colvar_set) @@ -321,6 +326,7 @@ CONTAINS NULLIFY(constraint_info%colv_intermolecular) NULLIFY(constraint_info%colv_restraint) NULLIFY(constraint_info%colv_k0) + NULLIFY(constraint_info%colv_exclude_qm, constraint_info%colv_exclude_mm) ! G3x3 constraint_info%nconst_g33 = 0 NULLIFY(constraint_info%const_g33_mol) @@ -334,6 +340,7 @@ CONTAINS NULLIFY(constraint_info%g33_intermolecular) NULLIFY(constraint_info%g33_restraint) NULLIFY(constraint_info%g33_k0) + NULLIFY(constraint_info%g33_exclude_qm, constraint_info%g33_exclude_mm) ! G4x6 constraint_info%nconst_g46 = 0 NULLIFY(constraint_info%const_g46_mol) @@ -351,6 +358,7 @@ CONTAINS NULLIFY(constraint_info%g46_intermolecular) NULLIFY(constraint_info%g46_restraint) NULLIFY(constraint_info%g46_k0) + NULLIFY(constraint_info%g46_exclude_qm, constraint_info%g46_exclude_mm) END SUBROUTINE init_constraint @@ -608,6 +616,14 @@ CONTAINS DEALLOCATE(constraint_info%fixed_mol_k0, stat=stat) CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) END IF + IF(ASSOCIATED(constraint_info%fixed_exclude_qm)) THEN + DEALLOCATE(constraint_info%fixed_exclude_qm, stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + END IF + IF(ASSOCIATED(constraint_info%fixed_exclude_mm)) THEN + DEALLOCATE(constraint_info%fixed_exclude_mm, stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + END IF ! Collective Constraint IF (ASSOCIATED(constraint_info%colvar_set)) THEN DO i = 1, SIZE(constraint_info%colvar_set) @@ -643,6 +659,14 @@ CONTAINS DEALLOCATE(constraint_info%colv_k0,STAT=stat) CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) END IF + IF(ASSOCIATED(constraint_info%colv_exclude_qm)) THEN + DEALLOCATE(constraint_info%colv_exclude_qm, stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + END IF + IF(ASSOCIATED(constraint_info%colv_exclude_mm)) THEN + DEALLOCATE(constraint_info%colv_exclude_mm, stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + END IF ! G3x3 IF(ASSOCIATED(constraint_info%const_g33_mol)) THEN DEALLOCATE(constraint_info%const_g33_mol,STAT=stat) @@ -688,6 +712,14 @@ CONTAINS DEALLOCATE(constraint_info%g33_k0,STAT=stat) CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) END IF + IF(ASSOCIATED(constraint_info%g33_exclude_qm)) THEN + DEALLOCATE(constraint_info%g33_exclude_qm, stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + END IF + IF(ASSOCIATED(constraint_info%g33_exclude_mm)) THEN + DEALLOCATE(constraint_info%g33_exclude_mm, stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + END IF ! G4x6 IF(ASSOCIATED(constraint_info%const_g46_mol)) THEN DEALLOCATE(constraint_info%const_g46_mol,STAT=stat) @@ -749,6 +781,14 @@ CONTAINS DEALLOCATE(constraint_info%g46_k0,STAT=stat) CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) END IF + IF(ASSOCIATED(constraint_info%g46_exclude_qm)) THEN + DEALLOCATE(constraint_info%g46_exclude_qm, stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + END IF + IF(ASSOCIATED(constraint_info%g46_exclude_mm)) THEN + DEALLOCATE(constraint_info%g46_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/QMMM/SE/regtest_2/TEST_FILES b/tests/QMMM/SE/regtest_2/TEST_FILES new file mode 100644 index 0000000000..86ea2a2ba8 --- /dev/null +++ b/tests/QMMM/SE/regtest_2/TEST_FILES @@ -0,0 +1,14 @@ +# runs are executed in the same order as in this file +# the second field tells which test should be run in order to compare with the last available output +# e.g. 0 means do not compare anything, running is enough +# 1 compares the last total energy in the file +# for details see cp2k/tools/do_regtest +# QM/MM +nh3_g4x6_excl_mm.inp 2 +nh3_g4x6_excl_qm.inp 2 +water_colv_excl_mm.inp 2 +water_colv_excl_qm.inp 2 +water_fixd_excl_mm.inp 2 +water_fixd_excl_qm.inp 2 +water_g3x3_excl_mm.inp 2 +water_g3x3_excl_qm.inp 2 diff --git a/tests/QMMM/SE/regtest_2/TEST_FILES_RESET b/tests/QMMM/SE/regtest_2/TEST_FILES_RESET new file mode 100644 index 0000000000..792d600548 --- /dev/null +++ b/tests/QMMM/SE/regtest_2/TEST_FILES_RESET @@ -0,0 +1 @@ +# diff --git a/tests/QMMM/SE/regtest_2/nh3_g4x6_excl_mm.inp b/tests/QMMM/SE/regtest_2/nh3_g4x6_excl_mm.inp new file mode 100644 index 0000000000..83ca94b215 --- /dev/null +++ b/tests/QMMM/SE/regtest_2/nh3_g4x6_excl_mm.inp @@ -0,0 +1,91 @@ +&FORCE_EVAL + METHOD QMMM + &DFT + &QS + METHOD AM1 + &END QS + &SCF + SCF_GUESS ATOMIC + &END SCF + &END DFT + &QMMM + NOCENTER + &CELL + ABC 6.0 6.0 6.0 + UNIT ANGSTROM + &END CELL + ECOUPL NONE + &QM_KIND H + MM_INDEX 2 3 4 + &END QM_KIND + &QM_KIND N + MM_INDEX 1 + &END QM_KIND + &END + &MM + &FORCEFIELD + parm_file_name ../../../Fist/sample_pot/nh3.pot + parmtype CHM + &END FORCEFIELD + &POISSON + &EWALD + EWALD_TYPE spme + ALPHA .44 + GMAX 30 + O_SPLINE 6 + &END EWALD + &END POISSON + &END MM + &SUBSYS + &CELL + ABC 10.0 10.0 10.0 + UNIT ANGSTROM + &END CELL + &TOPOLOGY + COORD_FILE_NAME ../sample_pdb/NH3.pdb + COORDINATE pdb + CONNECTIVITY MOL_SET + &MOL_SET + &MOLECULE + NMOL 2 + CONN_FILE_NAME ../../../Fist/sample_psf/nh3.psf + CONNECTIVITY UPSF + &END + &END + &END TOPOLOGY + &PRINT + &MOLECULES + &END + &END + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + FFTLIB FFTSG + PROJECT NH3-G46-EXMM + RUN_TYPE md + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CONSTRAINT + CONSTRAIN_INIT T + &G4X6 + DISTANCES 1.9275 1.9275 1.9275 3.1483 3.1483 3.1483 + MOLNAME MOL01 + ATOMS 1 2 3 4 + EXCLUDE_MM + &END G4X6 + &END CONSTRAINT + &MD + ENSEMBLE NVE + STEPS 20 + TIMESTEP 0.5 + TEMPERATURE 298 + &NOSE + LENGTH 3 + YOSHIDA 3 + TIMECON 1000 + MTS 2 + &END NOSE + &END MD +&END MOTION + diff --git a/tests/QMMM/SE/regtest_2/nh3_g4x6_excl_qm.inp b/tests/QMMM/SE/regtest_2/nh3_g4x6_excl_qm.inp new file mode 100644 index 0000000000..07663c7cc8 --- /dev/null +++ b/tests/QMMM/SE/regtest_2/nh3_g4x6_excl_qm.inp @@ -0,0 +1,91 @@ +&FORCE_EVAL + METHOD QMMM + &DFT + &QS + METHOD AM1 + &END QS + &SCF + SCF_GUESS ATOMIC + &END SCF + &END DFT + &QMMM + NOCENTER + &CELL + ABC 6.0 6.0 6.0 + UNIT ANGSTROM + &END CELL + ECOUPL NONE + &QM_KIND H + MM_INDEX 2 3 4 + &END QM_KIND + &QM_KIND N + MM_INDEX 1 + &END QM_KIND + &END + &MM + &FORCEFIELD + parm_file_name ../../../Fist/sample_pot/nh3.pot + parmtype CHM + &END FORCEFIELD + &POISSON + &EWALD + EWALD_TYPE spme + ALPHA .44 + GMAX 30 + O_SPLINE 6 + &END EWALD + &END POISSON + &END MM + &SUBSYS + &CELL + ABC 10.0 10.0 10.0 + UNIT ANGSTROM + &END CELL + &TOPOLOGY + COORD_FILE_NAME ../sample_pdb/NH3.pdb + COORDINATE pdb + CONNECTIVITY MOL_SET + &MOL_SET + &MOLECULE + NMOL 2 + CONN_FILE_NAME ../../../Fist/sample_psf/nh3.psf + CONNECTIVITY UPSF + &END + &END + &END TOPOLOGY + &PRINT + &MOLECULES + &END + &END + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + FFTLIB FFTSG + PROJECT NH3-G46-EXQM + RUN_TYPE md + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CONSTRAINT + CONSTRAIN_INIT T + &G4X6 + DISTANCES 1.9275 1.9275 1.9275 3.1483 3.1483 3.1483 + MOLNAME MOL01 + ATOMS 1 2 3 4 + EXCLUDE_QM + &END G4X6 + &END CONSTRAINT + &MD + ENSEMBLE NVE + STEPS 20 + TIMESTEP 0.5 + TEMPERATURE 298 + &NOSE + LENGTH 3 + YOSHIDA 3 + TIMECON 1000 + MTS 2 + &END NOSE + &END MD +&END MOTION + diff --git a/tests/QMMM/SE/regtest_2/tidy b/tests/QMMM/SE/regtest_2/tidy new file mode 100644 index 0000000000..8e44cbdbc5 --- /dev/null +++ b/tests/QMMM/SE/regtest_2/tidy @@ -0,0 +1,18 @@ +#!/bin/tcsh +# Script to clean files after execution of a test +echo "This script will clean the directory" $PWD "from scratch files.." +echo +foreach file (`cvs update -A -P|grep '? '|awk '{print $2}'`) + if ( ${file:e} == 'inp' ) then + echo "This files could be a new input file you don't want to be deleted.. Do you really want it deleted?" + if ( $1 == "-int" ) then + echo "You requested to run the cleaning script interactively! Forget about complains!" + rm -f $file + else + rm -i $file + endif + else + rm -f $file + endif +end +rm -f *.bak diff --git a/tests/QMMM/SE/regtest_2/water_colv_excl_mm.inp b/tests/QMMM/SE/regtest_2/water_colv_excl_mm.inp new file mode 100644 index 0000000000..13e7146aff --- /dev/null +++ b/tests/QMMM/SE/regtest_2/water_colv_excl_mm.inp @@ -0,0 +1,93 @@ + +&FORCE_EVAL + METHOD QMMM + &DFT + &QS + METHOD AM1 + &END QS + &SCF + SCF_GUESS ATOMIC + &END SCF + &END DFT + &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 .44 + GMAX 30 + &END EWALD + &END POISSON + &END MM + &QMMM + &CELL + ABC 6.0 6.0 6.0 + UNIT ANGSTROM + &END CELL + ECOUPL NONE + &QM_KIND H + MM_INDEX 2 3 + &END QM_KIND + &QM_KIND O + MM_INDEX 1 + &END QM_KIND + &END + &SUBSYS + &CELL + ABC 24.955 24.955 24.955 + UNIT ANGSTROM + &END CELL + &TOPOLOGY + COORD_FILE_NAME ../../../Fist/sample_pdb/water_2.pdb + COORDINATE pdb + &END TOPOLOGY + &PRINT + &MOLECULES + &END + &END + &COLVAR + &DISTANCE + ATOMS 2 3 + &END + &END + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + FFTLIB FFTSG + PROJECT WAT-COLV-EXMM + RUN_TYPE MD +&END GLOBAL +&MOTION + &CONSTRAINT + CONSTRAIN_INIT T + &COLLECTIVE + COLVAR 1 + TARGET 3.0 + MOLNAME WAT + EXCLUDE_MM + &END + &END + &MD + ENSEMBLE NVE + STEPS 20 + TIMESTEP 0.5 + TEMPERATURE 298 + &NOSE + LENGTH 3 + YOSHIDA 3 + TIMECON 1000 + MTS 2 + &END NOSE + &END MD +&END MOTION diff --git a/tests/QMMM/SE/regtest_2/water_colv_excl_qm.inp b/tests/QMMM/SE/regtest_2/water_colv_excl_qm.inp new file mode 100644 index 0000000000..15de2b6962 --- /dev/null +++ b/tests/QMMM/SE/regtest_2/water_colv_excl_qm.inp @@ -0,0 +1,93 @@ + +&FORCE_EVAL + METHOD QMMM + &DFT + &QS + METHOD AM1 + &END QS + &SCF + SCF_GUESS ATOMIC + &END SCF + &END DFT + &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 .44 + GMAX 30 + &END EWALD + &END POISSON + &END MM + &QMMM + &CELL + ABC 6.0 6.0 6.0 + UNIT ANGSTROM + &END CELL + ECOUPL NONE + &QM_KIND H + MM_INDEX 2 3 + &END QM_KIND + &QM_KIND O + MM_INDEX 1 + &END QM_KIND + &END + &SUBSYS + &CELL + ABC 24.955 24.955 24.955 + UNIT ANGSTROM + &END CELL + &TOPOLOGY + COORD_FILE_NAME ../../../Fist/sample_pdb/water_2.pdb + COORDINATE pdb + &END TOPOLOGY + &PRINT + &MOLECULES + &END + &END + &COLVAR + &DISTANCE + ATOMS 2 3 + &END + &END + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + FFTLIB FFTSG + PROJECT WAT-COLV-EXQM + RUN_TYPE MD +&END GLOBAL +&MOTION + &CONSTRAINT + CONSTRAIN_INIT T + &COLLECTIVE + COLVAR 1 + TARGET 3.0 + MOLNAME WAT + EXCLUDE_QM + &END + &END + &MD + ENSEMBLE NVE + STEPS 20 + TIMESTEP 0.5 + TEMPERATURE 298 + &NOSE + LENGTH 3 + YOSHIDA 3 + TIMECON 1000 + MTS 2 + &END NOSE + &END MD +&END MOTION diff --git a/tests/QMMM/SE/regtest_2/water_fixd_excl_mm.inp b/tests/QMMM/SE/regtest_2/water_fixd_excl_mm.inp new file mode 100644 index 0000000000..a7e6e7580f --- /dev/null +++ b/tests/QMMM/SE/regtest_2/water_fixd_excl_mm.inp @@ -0,0 +1,87 @@ + +&FORCE_EVAL + METHOD QMMM + &DFT + &QS + METHOD AM1 + &END QS + &SCF + SCF_GUESS ATOMIC + &END SCF + &END DFT + &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 .44 + GMAX 30 + &END EWALD + &END POISSON + &END MM + &QMMM + NOCENTER + &CELL + ABC 6.0 6.0 6.0 + UNIT ANGSTROM + &END CELL + ECOUPL NONE + &QM_KIND H + MM_INDEX 2 3 + &END QM_KIND + &QM_KIND O + MM_INDEX 1 + &END QM_KIND + &END + &SUBSYS + &CELL + ABC 24.955 24.955 24.955 + UNIT ANGSTROM + &END CELL + &TOPOLOGY + COORD_FILE_NAME ../../../Fist/sample_pdb/water_2.pdb + COORDINATE pdb + &END TOPOLOGY + &PRINT + &MOLECULES + &END + &END + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + FFTLIB FFTSG + PROJECT WAT-FXD-EXMM + RUN_TYPE MD +&END GLOBAL +&MOTION + &CONSTRAINT + CONSTRAIN_INIT T + &FIXED_ATOMS + MOLNAME WAT + EXCLUDE_MM + &END + &END + &MD + ENSEMBLE NVE + STEPS 20 + TIMESTEP 0.5 + TEMPERATURE 298 + &NOSE + LENGTH 3 + YOSHIDA 3 + TIMECON 1000 + MTS 2 + &END NOSE + &END MD +&END MOTION diff --git a/tests/QMMM/SE/regtest_2/water_fixd_excl_qm.inp b/tests/QMMM/SE/regtest_2/water_fixd_excl_qm.inp new file mode 100644 index 0000000000..8c2ed9ca0e --- /dev/null +++ b/tests/QMMM/SE/regtest_2/water_fixd_excl_qm.inp @@ -0,0 +1,87 @@ + +&FORCE_EVAL + METHOD QMMM + &DFT + &QS + METHOD AM1 + &END QS + &SCF + SCF_GUESS ATOMIC + &END SCF + &END DFT + &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 .44 + GMAX 30 + &END EWALD + &END POISSON + &END MM + &QMMM + NOCENTER + &CELL + ABC 6.0 6.0 6.0 + UNIT ANGSTROM + &END CELL + ECOUPL NONE + &QM_KIND H + MM_INDEX 2 3 + &END QM_KIND + &QM_KIND O + MM_INDEX 1 + &END QM_KIND + &END + &SUBSYS + &CELL + ABC 24.955 24.955 24.955 + UNIT ANGSTROM + &END CELL + &TOPOLOGY + COORD_FILE_NAME ../../../Fist/sample_pdb/water_2.pdb + COORDINATE pdb + &END TOPOLOGY + &PRINT + &MOLECULES + &END + &END + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + FFTLIB FFTSG + PROJECT WAT-FXD-EXQM + RUN_TYPE MD +&END GLOBAL +&MOTION + &CONSTRAINT + CONSTRAIN_INIT T + &FIXED_ATOMS + MOLNAME WAT + EXCLUDE_QM + &END + &END + &MD + ENSEMBLE NVE + STEPS 20 + TIMESTEP 0.5 + TEMPERATURE 298 + &NOSE + LENGTH 3 + YOSHIDA 3 + TIMECON 1000 + MTS 2 + &END NOSE + &END MD +&END MOTION diff --git a/tests/QMMM/SE/regtest_2/water_g3x3_excl_mm.inp b/tests/QMMM/SE/regtest_2/water_g3x3_excl_mm.inp new file mode 100644 index 0000000000..6dcb61838a --- /dev/null +++ b/tests/QMMM/SE/regtest_2/water_g3x3_excl_mm.inp @@ -0,0 +1,88 @@ + +&FORCE_EVAL + METHOD QMMM + &DFT + &QS + METHOD AM1 + &END QS + &SCF + SCF_GUESS ATOMIC + &END SCF + &END DFT + &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 .44 + GMAX 30 + &END EWALD + &END POISSON + &END MM + &QMMM + &CELL + ABC 6.0 6.0 6.0 + UNIT ANGSTROM + &END CELL + ECOUPL NONE + &QM_KIND H + MM_INDEX 2 3 + &END QM_KIND + &QM_KIND O + MM_INDEX 1 + &END QM_KIND + &END + &SUBSYS + &CELL + ABC 24.955 24.955 24.955 + UNIT ANGSTROM + &END CELL + &TOPOLOGY + COORD_FILE_NAME ../../../Fist/sample_pdb/water_2.pdb + COORDINATE pdb + &END TOPOLOGY + &PRINT + &MOLECULES + &END + &END + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + FFTLIB FFTSG + PROJECT WAT-G33-EXMM + RUN_TYPE MD +&END GLOBAL +&MOTION + &CONSTRAINT + CONSTRAIN_INIT T + &G3X3 + DISTANCES 1.8897268 1.8897268 3.0859239 + MOLNAME WAT + ATOMS 1 2 3 + EXCLUDE_MM + &END G3X3 + &END + &MD + ENSEMBLE NVE + STEPS 20 + TIMESTEP 0.5 + TEMPERATURE 298 + &NOSE + LENGTH 3 + YOSHIDA 3 + TIMECON 1000 + MTS 2 + &END NOSE + &END MD +&END MOTION diff --git a/tests/QMMM/SE/regtest_2/water_g3x3_excl_qm.inp b/tests/QMMM/SE/regtest_2/water_g3x3_excl_qm.inp new file mode 100644 index 0000000000..689eb3aff0 --- /dev/null +++ b/tests/QMMM/SE/regtest_2/water_g3x3_excl_qm.inp @@ -0,0 +1,88 @@ + +&FORCE_EVAL + METHOD QMMM + &DFT + &QS + METHOD AM1 + &END QS + &SCF + SCF_GUESS ATOMIC + &END SCF + &END DFT + &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 .44 + GMAX 30 + &END EWALD + &END POISSON + &END MM + &QMMM + &CELL + ABC 6.0 6.0 6.0 + UNIT ANGSTROM + &END CELL + ECOUPL NONE + &QM_KIND H + MM_INDEX 2 3 + &END QM_KIND + &QM_KIND O + MM_INDEX 1 + &END QM_KIND + &END + &SUBSYS + &CELL + ABC 24.955 24.955 24.955 + UNIT ANGSTROM + &END CELL + &TOPOLOGY + COORD_FILE_NAME ../../../Fist/sample_pdb/water_2.pdb + COORDINATE pdb + &END TOPOLOGY + &PRINT + &MOLECULES + &END + &END + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + FFTLIB FFTSG + PROJECT WAT-G33-EXQM + RUN_TYPE MD +&END GLOBAL +&MOTION + &CONSTRAINT + CONSTRAIN_INIT T + &G3X3 + DISTANCES 1.8897268 1.8897268 3.0859239 + MOLNAME WAT + ATOMS 1 2 3 + EXCLUDE_QM + &END G3X3 + &END + &MD + ENSEMBLE NVE + STEPS 20 + TIMESTEP 0.5 + TEMPERATURE 298 + &NOSE + LENGTH 3 + YOSHIDA 3 + TIMECON 1000 + MTS 2 + &END NOSE + &END MD +&END MOTION diff --git a/tests/QMMM/SE/sample_pdb/NH3.pdb b/tests/QMMM/SE/sample_pdb/NH3.pdb new file mode 100644 index 0000000000..59591103e9 --- /dev/null +++ b/tests/QMMM/SE/sample_pdb/NH3.pdb @@ -0,0 +1,8 @@ +ATOM 5 NH3 NH3 A 2 -5.612 0.338 1.176 1.00 -1.14 AMMO N +ATOM 6 HC NH3 A 2 -4.751 0.181 0.654 1.00 0.38 AMMO H +ATOM 7 HC NH3 A 2 -5.594 -0.259 2.043 1.00 0.38 AMMO H +ATOM 8 HC NH3 A 2 -5.450 1.320 1.518 1.00 0.38 AMMO H +ATOM 21 NH3 NH3 A 6 -5.431 -5.410 1.645 1.00 -1.14 AMMO N +ATOM 22 HC NH3 A 6 -4.538 -5.350 1.105 1.00 0.38 AMMO H +ATOM 23 HC NH3 A 6 -5.088 -6.094 2.380 1.00 0.38 AMMO H +ATOM 24 HC NH3 A 6 -6.119 -5.889 1.072 1.00 0.38 AMMO H diff --git a/tests/TEST_DIRS b/tests/TEST_DIRS index c8df2be8c9..e8b12cc3bc 100644 --- a/tests/TEST_DIRS +++ b/tests/TEST_DIRS @@ -23,6 +23,7 @@ QS/regtest-gapw NEB/regtest-2 QS/regtest-nmr-uks-1 QMMM/SE/regtest +QMMM/SE/regtest_2 QS/regtest-gapw-2 Fist/regtest-4 QS/regtest-nmr-1