From ed13711e52bcde5115587b452195acbeb4fe0644 Mon Sep 17 00:00:00 2001 From: Teodoro Laino Date: Mon, 23 Nov 2009 14:22:14 +0000 Subject: [PATCH] out of plane bends + regtests (Louis) svn-origin-rev: 9215 --- src/fist_force.F | 76 ++--- src/fist_intra_force.F | 143 +++++++--- src/force_field_types.F | 99 ++++++- src/force_fields_all.F | 315 +++++++++++++++++++-- src/force_fields_input.F | 60 +++- src/force_fields_util.F | 177 +++++++++--- src/input_constants.F | 170 +++++------ src/input_cp2k_mm.F | 75 ++++- src/mol_force.F | 137 ++++++++- src/molecule_kind_types.F | 104 +++++-- src/topology_connectivity_util.F | 80 +++++- tests/Fist/regtest-14/TEST_FILES | 1 + tests/Fist/regtest-14/imp_test_12.inp | 5 +- tests/Fist/regtest-14/opbend_test.inp | 178 ++++++++++++ tests/Fist/sample_psf/bromochloroamine.psf | 28 ++ tests/Fist/sample_psf/imp_test.psf | 6 +- 16 files changed, 1377 insertions(+), 277 deletions(-) create mode 100644 tests/Fist/regtest-14/opbend_test.inp create mode 100644 tests/Fist/sample_psf/bromochloroamine.psf diff --git a/src/fist_force.F b/src/fist_force.F index 20411c6685..b3505e87e4 100644 --- a/src/fist_force.F +++ b/src/fist_force.F @@ -82,10 +82,13 @@ MODULE fist_force REAL (KIND=dp) :: pot_manybody, pot_bend, pot_torsion REAL (KIND=dp) :: pot_nonbond, pot_g, pot_bond REAL (KIND=dp) :: pot_imptors, pot_onef, pot_urey_bradley + REAL (KIND=dp) :: pot_opbend REAL (KIND=dp), DIMENSION (:,:), POINTER :: f_nonbond, f_g, f_bond, f_bend,& - f_torsion, f_imptors, f_onef, f_ub + f_torsion, f_imptors, f_onef, f_ub,& + f_opbend REAL (KIND=dp), DIMENSION (3,3) :: pv_nonbond, pv_g, pv_bond, pv_onef, pv_ub,& - pv_bend, pv_torsion, pv_imptors + pv_bend, pv_torsion, pv_imptors,& + pv_opbend END TYPE debug_variables_type CONTAINS @@ -123,13 +126,13 @@ CONTAINS failure, shell_model_ad, & shell_present, use_virial REAL(KIND=dp) :: mass, massc, masss, pot_bend, pot_bond, pot_imptors, & - pot_manybody, pot_nonbond, pot_onef, pot_shell, pot_torsion, & - pot_urey_bradley, umass, vg_coulomb, xdum1, xdum2, xdum3 + pot_manybody, pot_nonbond, pot_onef, pot_opbend, pot_shell, & + pot_torsion, pot_urey_bradley, umass, vg_coulomb, xdum1, xdum2, xdum3 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: f_nonbond, f_total, & fcore_nonbond, fcore_total, fg_coulomb, fgcore_coulomb, & fgshell_coulomb, fshell_nonbond, fshell_total REAL(KIND=dp), DIMENSION(3, 3) :: ident, pv_bc, pv_bend, pv_bond, pv_g, & - pv_imptors, pv_nonbond, pv_onef, pv_torsion, pv_urey_bradley + pv_imptors, pv_nonbond, pv_onef, pv_opbend, pv_torsion, pv_urey_bradley TYPE(atomic_kind_type), DIMENSION(:), & POINTER :: atomic_kind_set TYPE(atomic_kind_type), POINTER :: atomic_kind @@ -234,13 +237,14 @@ CONTAINS core_particle_set ( i ) % f_embed = 0.0_dp END DO ENDIF - + IF (use_virial) THEN pv_bc = 0.0_dp pv_bond = 0.0_dp pv_bend = 0.0_dp pv_torsion = 0.0_dp pv_imptors = 0.0_dp + pv_opbend = 0.0_dp pv_onef = 0.0_dp pv_urey_bradley = 0.0_dp pv_nonbond = 0.0_dp @@ -253,6 +257,7 @@ CONTAINS pot_bond = 0.0_dp pot_bend = 0.0_dp pot_torsion = 0.0_dp + pot_opbend = 0.0_dp pot_imptors = 0.0_dp pot_onef = 0.0_dp pot_urey_bradley = 0.0_dp @@ -267,7 +272,7 @@ CONTAINS WRITE(iw,'(3f15.9)')((particle_set ( i ) % f ( j ),j=1,3), i = 1, SIZE(particle_set)) END IF - IF (fist_nonbond_env%do_nonbonded) THEN + IF (fist_nonbond_env%do_nonbonded) THEN ! Compute density for EAM CALL density_nonbond ( fist_nonbond_env, particle_set, cell, para_env ,error=error) @@ -326,8 +331,8 @@ CONTAINS CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) ALLOCATE ( fgcore_coulomb(3, nshell), STAT=stat ) CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) - fgshell_coulomb = 0.0_dp - fgcore_coulomb = 0.0_dp + fgshell_coulomb = 0.0_dp + fgcore_coulomb = 0.0_dp END IF IF(shell_present.AND.do_multipoles) THEN CALL cp_unimplemented_error(fromWhere=routineP, & @@ -439,7 +444,7 @@ CONTAINS CALL bonded_correct_gaussian ( atomic_kind_set, local_particles, particle_set,& ewald_env, thermo%e_bonded, pv_bc, shell_particle_set=shell_particle_set,& core_particle_set=core_particle_set, cell=cell, use_virial=use_virial, error=error) - + ELSE IF (.NOT. do_multipoles) THEN CALL bonded_correct_gaussian ( atomic_kind_set, local_particles,& @@ -474,22 +479,23 @@ CONTAINS IF ( PRESENT ( debug ) ) THEN CALL force_intra_control ( molecule_set, molecule_kind_set, local_molecules, & particle_set, shell_particle_set, core_particle_set, pot_bond, pot_bend,& - pot_urey_bradley, pot_torsion, pot_imptors, pot_onef, pot_shell, pv_bond,& - pv_bend, pv_urey_bradley, pv_torsion, pv_imptors, pv_onef, debug%f_bond,& - debug%f_bend, debug%f_torsion, debug%f_ub, debug%f_imptors,debug%f_onef,& - cell, use_virial, error ) + pot_urey_bradley, pot_torsion, pot_imptors, pot_opbend, pot_onef, & + pot_shell, pv_bond, pv_bend, pv_urey_bradley, pv_torsion, pv_imptors, & + pv_opbend, pv_onef, debug%f_bond, debug%f_bend, debug%f_torsion, debug%f_ub, & + debug%f_imptors, debug%f_opbend, debug%f_onef, cell, use_virial, error ) + ELSE CALL force_intra_control ( molecule_set, molecule_kind_set, local_molecules, & particle_set, shell_particle_set, core_particle_set, pot_bond, pot_bend,& - pot_urey_bradley, pot_torsion, pot_imptors, pot_onef, pot_shell, pv_bond,& - pv_bend, pv_urey_bradley, pv_torsion, pv_imptors, pv_onef, cell=cell,& - use_virial=use_virial, error=error) + pot_urey_bradley, pot_torsion, pot_imptors, pot_opbend, pot_onef, & + pot_shell, pv_bond, pv_bend, pv_urey_bradley, pv_torsion, pv_imptors, & + pv_opbend, pv_onef, cell=cell, use_virial=use_virial, error=error) ENDIF IF (iw>0) THEN xdum1 = cp_unit_from_cp2k(pot_bond,"kcalmol",error=error) xdum2 = cp_unit_from_cp2k(pot_bend,"kcalmol",error=error) - xdum3 = cp_unit_from_cp2k(pot_torsion+pot_imptors,"kcalmol",error=error) + xdum3 = cp_unit_from_cp2k(pot_torsion+pot_imptors+pot_opbend,"kcalmol",error=error) WRITE(iw,'(1x,"BOND = ",f13.4,2x,"ANGLE = ",f13.4,2x,'//& '"DIHED = ",f13.4)')xdum1, xdum2, xdum3 xdum1 = cp_unit_from_cp2k(pot_onef,"kcalmol",error=error) @@ -509,8 +515,8 @@ CONTAINS END IF ! add up all the potential energies - thermo%pot = pot_nonbond + pot_bond + pot_bend + pot_torsion + pot_onef + & - pot_imptors + pot_urey_bradley + pot_manybody + pot_shell + thermo%pot = pot_nonbond + pot_bond + pot_bend + pot_torsion +pot_opbend + & + pot_onef + pot_imptors + pot_urey_bradley + pot_manybody + pot_shell CALL mp_sum ( thermo%pot, para_env%group ) IF(shell_present) THEN thermo%harm_shell = pot_shell @@ -598,7 +604,7 @@ CONTAINS WRITE(iw,'(3f15.9)')((fcore_total ( j, i ),j=1,3), i = 1, nshell) END IF END IF - + IF (use_virial) THEN ! Add up all the pressure tensors IF ( ewald_type == do_ewald_none ) THEN @@ -609,10 +615,10 @@ CONTAINS DO i = 1, 3 ident(i, i) = 1.0_dp END DO - + virial%pv_virial = pv_nonbond + pv_bond + pv_bend + pv_torsion + pv_imptors + pv_urey_bradley + pv_bc CALL mp_sum (virial%pv_virial, para_env%group) - + virial%pv_virial = virial%pv_virial + ident * thermo%e_neut virial%pv_virial = virial%pv_virial + pv_g END IF @@ -689,16 +695,16 @@ CONTAINS ELSE particle_set ( i )%f ( 1 ) = f_total ( 1, i ) particle_set ( i )%f ( 2 ) = f_total ( 2, i ) - particle_set ( i )%f ( 3 ) = f_total ( 3, i ) + particle_set ( i )%f ( 3 ) = f_total ( 3, i ) END IF END DO DO i = 1,nshell - shell_particle_set(i)%f(1) = fshell_total(1,i) - shell_particle_set(i)%f(2) = fshell_total(2,i) - shell_particle_set(i)%f(3) = fshell_total(3,i) - core_particle_set(i)%f(1) = fcore_total(1,i) - core_particle_set(i)%f(2) = fcore_total(2,i) - core_particle_set(i)%f(3) = fcore_total(3,i) + shell_particle_set(i)%f(1) = fshell_total(1,i) + shell_particle_set(i)%f(2) = fshell_total(2,i) + shell_particle_set(i)%f(3) = fshell_total(3,i) + core_particle_set(i)%f(1) = fcore_total(1,i) + core_particle_set(i)%f(2) = fcore_total(2,i) + core_particle_set(i)%f(3) = fcore_total(3,i) END DO ELSEIF(shell_present .AND. .NOT. shell_model_ad) THEN CALL cp_unimplemented_error(fromWhere=routineP, & @@ -740,6 +746,8 @@ CONTAINS debug%pot_torsion = pot_torsion CALL mp_sum ( pot_imptors, para_env%group ) debug%pot_imptors = pot_imptors + CALL mp_sum ( pot_opbend, para_env%group ) + debug%pot_opbend = pot_opbend CALL mp_sum ( pot_onef, para_env%group ) debug%pot_onef = pot_onef CALL mp_sum ( pot_urey_bradley, para_env%group ) @@ -844,7 +852,7 @@ END SUBROUTINE fist_force_control !> possibly using the berry phase formalism !> \par History !> [01.2006] created -!> [12.2007] tlaino - University of Zurich - debug and extended +!> [12.2007] tlaino - University of Zurich - debug and extended !> \author Teodoro Laino ! ***************************************************************************** SUBROUTINE fist_dipole(fist_env, print_section, atomic_kind_set, particle_set, cell, unit_nr, error) @@ -880,7 +888,7 @@ END SUBROUTINE fist_force_control description='[DIPOLE]' CALL section_vals_val_get(print_section,"DIPOLE%REF_POINT",r_vals=ref_point,error=error) CALL section_vals_val_get(print_section,"DIPOLE%PERIODIC",l_val=do_berry,error=error) - + CALL get_reference_point(rcc,drcc,fist_env=fist_env,reference=reference,ref_point=ref_point,error=error) ! Dipole deriv will be the derivative of the Dipole (dM/dt=\sum e_j v_j) @@ -896,7 +904,7 @@ END SUBROUTINE fist_force_control END DO ria = twopi * MATMUL ( cell%h_inv, rcc ) zphase = CMPLX ( COS(ria), SIN(ria), dp )**charge_tot - + dria = twopi * MATMUL ( cell%h_inv, drcc ) dzphase = charge_tot * CMPLX ( -SIN(ria), COS(ria), dp )**(charge_tot-1.0_dp) * dria @@ -905,7 +913,7 @@ END SUBROUTINE fist_force_control DO iparticle_kind = 1, SIZE(atomic_kind_set) atomic_kind => atomic_kind_set(iparticle_kind) CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=q, atom_list=atom_list) - + DO i = 1, SIZE(atom_list) ria = particle_set(atom_list(i))%r(:) ria = pbc(ria,cell) diff --git a/src/fist_intra_force.F b/src/fist_intra_force.F index 0f8694d8b1..a77283c0e0 100644 --- a/src/fist_intra_force.F +++ b/src/fist_intra_force.F @@ -20,13 +20,15 @@ MODULE fist_intra_force USE mol_force, ONLY: force_bends,& force_bonds,& force_imp_torsions,& + force_opbends,& force_torsions,& get_pv_bend,& get_pv_bond,& get_pv_torsion USE molecule_kind_types, ONLY: & bend_type, bond_type, get_molecule_kind, impr_type, & - molecule_kind_type, onfo_type, shell_type, torsion_type, ub_type + molecule_kind_type, onfo_type, opbend_type, shell_type, torsion_type, & + ub_type USE molecule_types_new, ONLY: get_molecule,& molecule_type USE particle_types, ONLY: particle_type @@ -53,9 +55,10 @@ CONTAINS ! ***************************************************************************** SUBROUTINE force_intra_control ( molecule_set, molecule_kind_set, local_molecules,& particle_set, shell_particle_set, core_particle_set, pot_bond, pot_bend,& - pot_urey_bradley, pot_torsion, pot_imp_torsion, pot_onef, pot_shell, pv_bond,& - pv_bend, pv_urey_bradley, pv_torsion, pv_imp_torsion, pv_onef, f_bond,& - f_bend, f_torsion, f_ub, f_imptor, f_onef, cell, use_virial, error ) + pot_urey_bradley, pot_torsion, pot_imp_torsion, pot_opbend, pot_onef, & + pot_shell, pv_bond, pv_bend, pv_urey_bradley, pv_torsion, pv_imp_torsion, & + pv_opbend, pv_onef, f_bond, f_bend, f_torsion, f_ub, f_imptor, f_opbend, & + f_onef, cell, use_virial, error ) TYPE(molecule_type), POINTER :: molecule_set( : ) TYPE(molecule_kind_type), POINTER :: molecule_kind_set( : ) @@ -63,17 +66,17 @@ CONTAINS TYPE(particle_type), POINTER :: particle_set( : ), & shell_particle_set( : ), & core_particle_set( : ) - REAL(KIND=dp), INTENT(INOUT) :: pot_bond, pot_bend, & - pot_urey_bradley, & - pot_torsion, pot_imp_torsion, & - pot_onef, pot_shell + REAL(KIND=dp), INTENT(INOUT) :: pot_bond, pot_bend, pot_urey_bradley, & + pot_torsion, pot_imp_torsion, pot_opbend, pot_onef, pot_shell REAL(KIND=dp), DIMENSION(:, :), & INTENT(INOUT) :: pv_bond, pv_bend, & pv_urey_bradley, pv_torsion, & - pv_imp_torsion, pv_onef + pv_imp_torsion, pv_opbend, & + pv_onef REAL(KIND=dp), DIMENSION(:, :), & INTENT(INOUT), OPTIONAL :: f_bond, f_bend, f_torsion, & - f_ub, f_imptor, f_onef + f_ub, f_imptor, f_opbend, & + f_onef TYPE(cell_type), POINTER :: cell LOGICAL, INTENT(IN) :: use_virial TYPE(cp_error_type), INTENT(inout) :: error @@ -82,14 +85,15 @@ CONTAINS routineP = moduleN//':'//routineN INTEGER :: first_atom, handle, i, ibend, ibond, ikind, imol, imul, & - index_a, index_b, index_c, index_d, ionef, ishell, itorsion, nbends, & - nbonds, nimptors, nkind, nmol_per_kind, nonefs, nshell, ntorsions, nub + index_a, index_b, index_c, index_d, ionef, iopbend, ishell, itorsion, & + nbends, nbonds, nimptors, nkind, nmol_per_kind, nonefs, nopbends, & + nshell, ntorsions, nub REAL(KIND=dp) :: d12, d32, dist, dist1, dist2, energy, fscalar, id12, & id32, is32, ism, isn, r0, rab2, rinv, rinv2, s32, sm, sn, theta REAL(KIND=dp), DIMENSION(3) :: b12, b32, g1, g2, g3, gt1, & gt2, gt3, gt4, k1, k2, k3, & - k4, rij, t12, t32, t34, t43, & - tm, tn + k4, rij, t12, t32, t34, t41, & + t42, t43, tm, tn TYPE(bend_type), POINTER :: bend_list( : ) TYPE(bond_type), POINTER :: bond_list( : ) TYPE(cp_logger_type), POINTER :: logger @@ -97,6 +101,7 @@ CONTAINS TYPE(molecule_kind_type), POINTER :: molecule_kind TYPE(molecule_type), POINTER :: molecule TYPE(onfo_type), POINTER :: onef_list( : ) + TYPE(opbend_type), POINTER :: opbend_list( : ) TYPE(shell_type), POINTER :: shell_list( : ) TYPE(spline_data_p_type), DIMENSION(:), & POINTER :: spl_p @@ -114,15 +119,16 @@ CONTAINS IF ( PRESENT ( f_imptor)) f_imptor = 0.0_dp IF ( PRESENT ( f_ub)) f_ub = 0.0_dp IF ( PRESENT ( f_onef)) f_onef = 0.0_dp - + pot_bond = 0.0_dp pot_bend = 0.0_dp pot_urey_bradley = 0.0_dp pot_torsion = 0.0_dp pot_imp_torsion = 0.0_dp + pot_opbend =0.0_dp pot_onef = 0.0_dp pot_shell = 0.0_dp - + nkind = SIZE ( molecule_kind_set ) MOL: DO ikind = 1, nkind nmol_per_kind = local_molecules % n_el ( ikind ) @@ -133,11 +139,11 @@ CONTAINS molecule_kind => molecule % molecule_kind CALL get_molecule_kind ( molecule_kind, nbend = nbends, nbond = nbonds, & nimpr = nimptors, nub = nub, ntorsion = ntorsions, & - nonfo = nonefs, nshell=nshell, & + nopbend = nopbends, nonfo = nonefs, nshell=nshell, & bond_list = bond_list, ub_list = ub_list, & bend_list = bend_list, torsion_list = torsion_list, & - impr_list = impr_list, onfo_list = onef_list,& - shell_list=shell_list ) + impr_list = impr_list, opbend_list = opbend_list, & + onfo_list = onef_list, shell_list=shell_list ) CALL get_molecule ( molecule, first_atom=first_atom ) @@ -152,7 +158,7 @@ CONTAINS bond_list ( ibond ) % bond_kind % cs, & energy, fscalar) pot_bond = pot_bond + energy - + particle_set ( index_a )% f(1) = particle_set ( index_a ) % f(1) -rij(1) * fscalar particle_set ( index_a )% f(2) = particle_set ( index_a ) % f(2) -rij(2) * fscalar particle_set ( index_a )% f(3) = particle_set ( index_a ) % f(3) -rij(3) * fscalar @@ -173,7 +179,7 @@ CONTAINS f_bond(2,index_b) = f_bond(2,index_b) + rij(2)*fscalar f_bond(3,index_b) = f_bond(3,index_b) + rij(3)*fscalar END IF - + END DO BOND SHELL: DO ishell = 1, nshell @@ -185,7 +191,7 @@ CONTAINS k1(1) = shell_list ( ishell ) % shell_kind % k_spring CALL force_bonds(do_ff_harmonic, rij, 0.0_dp, k1, 0.0_dp, energy, fscalar) pot_shell = pot_shell + energy - + core_particle_set ( index_b )% f(1) = core_particle_set ( index_b ) % f(1) -rij(1) * fscalar core_particle_set ( index_b )% f(2) = core_particle_set ( index_b ) % f(2) -rij(2) * fscalar core_particle_set ( index_b )% f(3) = core_particle_set ( index_b ) % f(3) -rij(3) * fscalar @@ -220,9 +226,9 @@ CONTAINS f_ub(:,index_a) = f_ub(:,index_a) - rij*fscalar f_ub(:,index_b) = f_ub(:,index_b) + rij*fscalar END IF - + END DO UREY_BRADLEY - + BEND: DO ibend = 1, nbends index_a = bend_list ( ibend ) % a + first_atom - 1 index_b = bend_list ( ibend ) % b + first_atom - 1 @@ -261,7 +267,7 @@ CONTAINS particle_set ( index_c )% f(1) = particle_set ( index_c ) % f(1) + g3(1) * fscalar particle_set ( index_c )% f(2) = particle_set ( index_c ) % f(2) + g3(2) * fscalar particle_set ( index_c )% f(3) = particle_set ( index_c ) % f(3) + g3(3) * fscalar - + ! computing the pressure tensor k1 = fscalar*g1 k3 = fscalar*g3 @@ -279,7 +285,7 @@ CONTAINS index_a = torsion_list ( itorsion ) % a + first_atom - 1 index_b = torsion_list ( itorsion ) % b + first_atom - 1 index_c = torsion_list ( itorsion ) % c + first_atom - 1 - index_d = torsion_list ( itorsion ) % d + first_atom - 1 + index_d = torsion_list ( itorsion ) % d + first_atom - 1 t12 = particle_set ( index_a ) % r - particle_set ( index_b ) % r t32 = particle_set ( index_c ) % r - particle_set ( index_b ) % r t34 = particle_set ( index_c ) % r - particle_set ( index_d ) % r @@ -324,13 +330,13 @@ CONTAINS particle_set ( index_d )% f(1) = particle_set ( index_d ) % f(1) + gt4(1) * fscalar particle_set ( index_d )% f(2) = particle_set ( index_d ) % f(2) + gt4(2) * fscalar particle_set ( index_d )% f(3) = particle_set ( index_d ) % f(3) + gt4(3) * fscalar - + ! computing the pressure tensor k1 = fscalar*gt1 k3 = fscalar*gt3 k4 = fscalar*gt4 IF (use_virial) CALL get_pv_torsion(k1,k3,k4,t12,t32,t43,pv_torsion) - + ! the contribution from the torsions. ONLY FOR DEBUG IF ( PRESENT ( f_torsion)) THEN f_torsion(:,index_a) = f_torsion(:,index_a) + fscalar*gt1 @@ -340,7 +346,7 @@ CONTAINS END IF END DO END DO TORSION - + IMP_TORSION: DO itorsion = 1, nimptors index_a = impr_list ( itorsion ) % a + first_atom - 1 index_b = impr_list ( itorsion ) % b + first_atom - 1 @@ -388,7 +394,7 @@ CONTAINS particle_set ( index_d )% f(1) = particle_set ( index_d ) % f(1) + gt4(1) * fscalar particle_set ( index_d )% f(2) = particle_set ( index_d ) % f(2) + gt4(2) * fscalar particle_set ( index_d )% f(3) = particle_set ( index_d ) % f(3) + gt4(3) * fscalar - + ! computing the pressure tensor k1 = fscalar*gt1 k3 = fscalar*gt3 @@ -403,7 +409,72 @@ CONTAINS f_imptor(:,index_d) = f_imptor(:,index_d) + fscalar*gt4 END IF END DO IMP_TORSION - + + OPBEND: DO iopbend = 1, nopbends + index_a = opbend_list ( iopbend ) % a + first_atom - 1 + index_b = opbend_list ( iopbend ) % b + first_atom - 1 + index_c = opbend_list ( iopbend ) % c + first_atom - 1 + index_d = opbend_list ( iopbend ) % d + first_atom - 1 + + t12 = particle_set ( index_a ) % r - particle_set ( index_b ) % r + t32 = particle_set ( index_c ) % r - particle_set ( index_b ) % r + t34 = particle_set ( index_c ) % r - particle_set ( index_d ) % r + t43 = particle_set ( index_d ) % r - particle_set ( index_c ) % r + t41 = particle_set ( index_d ) % r - particle_set ( index_a ) % r + t42 = pbc(t41+t12,cell) + t12 = pbc(t12,cell) + t32 = pbc(t32,cell) + t41 = pbc(t41,cell) + t43 = pbc(t43,cell) + ! tm = t32 x t12 + tm(1) = t32(2)*t12(3)-t12(2)*t32(3) + tm(2) = -t32(1)*t12(3)+t12(1)*t32(3) + tm(3) = t32(1)*t12(2)-t12(1)*t32(2) + ! tn = - t41 x tm + tn(1) = -( t41(2)*tm(3)-tm(2)*t41(3)) + tn(2) = -(-t41(1)*tm(3)+tm(1)*t41(3)) + tn(3) = -( t41(1)*tm(2)-tm(1)*t41(2)) + sm = SQRT ( DOT_PRODUCT(tm,tm) ) + ism = 1.0_dp/sm + sn = SQRT ( DOT_PRODUCT(tn,tn) ) + isn = 1.0_dp/sn + s32 = SQRT ( DOT_PRODUCT(t32,t32) ) + is32 = 1.0_dp/s32 + CALL force_opbends( opbend_list(iopbend)%opbend_kind%id_type, & + s32, is32,ism,isn,tm,tn, t41, t42, t43, & + opbend_list(iopbend)%opbend_kind%k, & + opbend_list(iopbend)%opbend_kind%phi0, & + gt1,gt2,gt3,gt4,energy,fscalar) + pot_opbend = pot_opbend + energy + particle_set ( index_a )% f(1) = particle_set ( index_a ) % f(1) + gt1(1) * fscalar + particle_set ( index_a )% f(2) = particle_set ( index_a ) % f(2) + gt1(2) * fscalar + particle_set ( index_a )% f(3) = particle_set ( index_a ) % f(3) + gt1(3) * fscalar + particle_set ( index_b )% f(1) = particle_set ( index_b ) % f(1) + gt2(1) * fscalar + particle_set ( index_b )% f(2) = particle_set ( index_b ) % f(2) + gt2(2) * fscalar + particle_set ( index_b )% f(3) = particle_set ( index_b ) % f(3) + gt2(3) * fscalar + particle_set ( index_c )% f(1) = particle_set ( index_c ) % f(1) + gt3(1) * fscalar + particle_set ( index_c )% f(2) = particle_set ( index_c ) % f(2) + gt3(2) * fscalar + particle_set ( index_c )% f(3) = particle_set ( index_c ) % f(3) + gt3(3) * fscalar + particle_set ( index_d )% f(1) = particle_set ( index_d ) % f(1) + gt4(1) * fscalar + particle_set ( index_d )% f(2) = particle_set ( index_d ) % f(2) + gt4(2) * fscalar + particle_set ( index_d )% f(3) = particle_set ( index_d ) % f(3) + gt4(3) * fscalar + + ! computing the pressure tensor + k1 = fscalar*gt1 + k3 = fscalar*gt3 + k4 = fscalar*gt4 + + IF (use_virial) CALL get_pv_torsion(k1,k3,k4,t12,t32,t43,pv_opbend) + + ! the contribution from the opbends. ONLY FOR DEBUG + IF ( PRESENT ( f_opbend)) THEN + f_opbend(:,index_a) = f_opbend(:,index_a) + fscalar*gt1 + f_opbend(:,index_b) = f_opbend(:,index_b) + fscalar*gt2 + f_opbend(:,index_c) = f_opbend(:,index_c) + fscalar*gt3 + f_opbend(:,index_d) = f_opbend(:,index_d) + fscalar*gt4 + END IF + END DO OPBEND + ONEF: DO ionef = 1, nonefs index_a = onef_list ( ionef ) % a + first_atom - 1 index_b = onef_list ( ionef ) % b + first_atom - 1 @@ -428,11 +499,11 @@ CONTAINS particle_set ( index_b )% f(1) = particle_set ( index_b ) % f(1) +rij(1) * fscalar particle_set ( index_b )% f(2) = particle_set ( index_b ) % f(2) +rij(2) * fscalar particle_set ( index_b )% f(3) = particle_set ( index_b ) % f(3) +rij(3) * fscalar - + ! computing the pressure tensor k2 = rij * fscalar IF (use_virial) CALL get_pv_bond ( k2,rij,pv_onef ) - + ! the contribution from the bonds. ONLY FOR DEBUG IF ( PRESENT ( f_onef ) ) THEN f_onef(:,index_a) = f_onef(:,index_a) - rij*fscalar @@ -441,10 +512,10 @@ CONTAINS END DO ONEF END DO END DO MOL - + CALL timestop(handle) - + END SUBROUTINE force_intra_control - + END MODULE fist_intra_force diff --git a/src/force_field_types.F b/src/force_field_types.F index 16886ca05a..43341fd0a8 100644 --- a/src/force_field_types.F +++ b/src/force_field_types.F @@ -32,7 +32,7 @@ MODULE force_field_types PRIVATE CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'force_field_types' - + ! *** Define the derived structure types *** ! ***************************************************************************** @@ -73,6 +73,13 @@ MODULE force_field_types INTEGER :: kind_number END TYPE impr_kind_type +! ***************************************************************************** + TYPE opbend_kind_type + INTEGER :: id_type + REAL(KIND = dp) :: k,phi0 + INTEGER :: kind_number + END TYPE opbend_kind_type + ! ***************************************************************************** TYPE onfo_kind_type TYPE ( spline_factor_type ), POINTER :: spl_f @@ -125,6 +132,13 @@ MODULE force_field_types CHARACTER (LEN=default_string_length), POINTER :: impr_d(:) REAL(KIND=dp), POINTER :: impr_k(:) REAL(KIND=dp), POINTER :: impr_phi0(:) + INTEGER, POINTER :: opbend_kind(:) + CHARACTER (LEN=default_string_length), POINTER :: opbend_a(:) + CHARACTER (LEN=default_string_length), POINTER :: opbend_b(:) + CHARACTER (LEN=default_string_length), POINTER :: opbend_c(:) + CHARACTER (LEN=default_string_length), POINTER :: opbend_d(:) + REAL(KIND=dp), POINTER :: opbend_k(:) + REAL(KIND=dp), POINTER :: opbend_phi0(:) TYPE (pair_potential_p_type), POINTER :: nonbonded TYPE (pair_potential_p_type), POINTER :: nonbonded14 TYPE (shell_p_type), DIMENSION(:), POINTER :: shell_list @@ -241,6 +255,7 @@ MODULE force_field_types allocate_ub_kind_set,& allocate_torsion_kind_set,& allocate_impr_kind_set,& + allocate_opbend_kind_set,& allocate_onfo_kind_set,& init_ff_type,& deallocate_ff_type,& @@ -259,6 +274,7 @@ MODULE force_field_types onfo_kind_type,& onfo_kind_dealloc_ref,& torsion_kind_type,& + opbend_kind_type,& ub_kind_type,& ub_kind_dealloc_ref,& charmm_info_type,& @@ -425,6 +441,34 @@ CONTAINS END DO END SUBROUTINE allocate_impr_kind_set +! ***************************************************************************** +!> \brief Allocate and initialize a opbend kind set. +! ***************************************************************************** + SUBROUTINE allocate_opbend_kind_set(opbend_kind_set,nkind,error) + + TYPE(opbend_kind_type), DIMENSION(:), & + POINTER :: opbend_kind_set + INTEGER, INTENT(IN) :: nkind + TYPE(cp_error_type), INTENT(inout) :: error + + CHARACTER(LEN=*), PARAMETER :: routineN = 'allocate_opbend_kind_set', & + routineP = moduleN//':'//routineN + + INTEGER :: ikind, stat + LOGICAL :: failure + + failure = .FALSE. + NULLIFY(opbend_kind_set) + ALLOCATE (opbend_kind_set(nkind),STAT=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + DO ikind=1,nkind + opbend_kind_set(ikind)%id_type = do_ff_undef + opbend_kind_set(ikind)%k = 0.0_dp + opbend_kind_set(ikind)%phi0 = 0.0_dp + opbend_kind_set(ikind)%kind_number = ikind + END DO + END SUBROUTINE allocate_opbend_kind_set + ! ***************************************************************************** !> \brief Allocate and initialize a onfo kind set. ! ***************************************************************************** @@ -469,7 +513,7 @@ CONTAINS failure = .FALSE. !----------------------------------------------------------------------------- - ! 1. Initialize + ! 1. Initialize !----------------------------------------------------------------------------- NULLIFY(ff_type%inp_info,ff_type%chm_info,ff_type%amb_info,ff_type%gro_info) ALLOCATE(ff_type%inp_info,STAT=stat) @@ -485,7 +529,7 @@ CONTAINS ! 2. Initialize and Nullify things in ff_type%inp_info !----------------------------------------------------------------------------- CALL init_inp_info(ff_type%inp_info) - + !----------------------------------------------------------------------------- ! 3. Initialize and Nullify things in ff_type%chm_info !----------------------------------------------------------------------------- @@ -500,7 +544,7 @@ CONTAINS ! 5. Initialize and Nullify things in ff_type%amb_info !----------------------------------------------------------------------------- CALL init_amber_info(ff_type%amb_info) - + END SUBROUTINE init_ff_type ! ***************************************************************************** @@ -658,6 +702,13 @@ CONTAINS NULLIFY(inp_info%impr_d) NULLIFY(inp_info%impr_k) NULLIFY(inp_info%impr_phi0) + NULLIFY(inp_info%opbend_kind) + NULLIFY(inp_info%opbend_a) + NULLIFY(inp_info%opbend_b) + NULLIFY(inp_info%opbend_c) + NULLIFY(inp_info%opbend_d) + NULLIFY(inp_info%opbend_k) + NULLIFY(inp_info%opbend_phi0) NULLIFY(inp_info%nonbonded) NULLIFY(inp_info%nonbonded14) NULLIFY(inp_info%shell_list) @@ -708,15 +759,15 @@ CONTAINS END IF IF(ASSOCIATED(ff_type%chm_info)) THEN DEALLOCATE(ff_type%chm_info,STAT=stat) - CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) END IF IF(ASSOCIATED(ff_type%gro_info)) THEN DEALLOCATE(ff_type%gro_info,STAT=stat) - CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) END IF IF(ASSOCIATED(ff_type%amb_info)) THEN DEALLOCATE(ff_type%amb_info,STAT=stat) - CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) END IF END SUBROUTINE deallocate_ff_type @@ -1244,6 +1295,34 @@ CONTAINS DEALLOCATE(inp_info%impr_phi0,STAT=stat) CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) END IF + IF(ASSOCIATED(inp_info%opbend_kind)) THEN + DEALLOCATE(inp_info%opbend_kind,STAT=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + END IF + IF(ASSOCIATED(inp_info%opbend_a)) THEN + DEALLOCATE(inp_info%opbend_a,STAT=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + END IF + IF(ASSOCIATED(inp_info%opbend_b)) THEN + DEALLOCATE(inp_info%opbend_b,STAT=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + END IF + IF(ASSOCIATED(inp_info%opbend_c)) THEN + DEALLOCATE(inp_info%opbend_c,STAT=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + END IF + IF(ASSOCIATED(inp_info%opbend_d)) THEN + DEALLOCATE(inp_info%opbend_d,STAT=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + END IF + IF(ASSOCIATED(inp_info%opbend_k)) THEN + DEALLOCATE(inp_info%opbend_k,STAT=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + END IF + IF(ASSOCIATED(inp_info%opbend_phi0)) THEN + DEALLOCATE(inp_info%opbend_phi0,STAT=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + END IF IF(ASSOCIATED(inp_info%nonbonded)) THEN CALL pair_potential_p_release(inp_info%nonbonded,error=error) END IF @@ -1325,7 +1404,7 @@ CONTAINS DEALLOCATE(torsion_kind%phi0,STAT=stat) CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) END IF - + END SUBROUTINE torsion_kind_dealloc_ref ! ***************************************************************************** @@ -1362,11 +1441,9 @@ CONTAINS LOGICAL :: failure failure = .FALSE. - ! ! Questa e' la migliore routine che mente umana abbia concepito! ;-) - ! - + ! END SUBROUTINE impr_kind_dealloc_ref ! ***************************************************************************** diff --git a/src/force_fields_all.F b/src/force_fields_all.F index e94e757edf..9e095db02d 100644 --- a/src/force_fields_all.F +++ b/src/force_fields_all.F @@ -26,10 +26,11 @@ MODULE force_fields_all USE force_field_types, ONLY: & allocate_bend_kind_set, allocate_bond_kind_set, & allocate_impr_kind_set, allocate_onfo_kind_set, & - allocate_torsion_kind_set, allocate_ub_kind_set, amber_info_type, & - bend_kind_type, bond_kind_type, charmm_info_type, force_field_type, & - gromos_info_type, impr_kind_type, input_info_type, onfo_kind_type, & - torsion_kind_type, ub_kind_type + allocate_opbend_kind_set, allocate_torsion_kind_set, & + allocate_ub_kind_set, amber_info_type, bend_kind_type, bond_kind_type, & + charmm_info_type, force_field_type, gromos_info_type, impr_kind_type, & + input_info_type, onfo_kind_type, opbend_kind_type, torsion_kind_type, & + ub_kind_type USE input_constants, ONLY: do_ff_amber,& do_ff_charmm,& do_ff_g87,& @@ -45,8 +46,8 @@ MODULE force_fields_all USE memory_utilities, ONLY: reallocate USE molecule_kind_types, ONLY: & bend_type, bond_type, get_molecule_kind, impr_type, & - molecule_kind_type, onfo_type, set_molecule_kind, shell_type, & - torsion_type, ub_type + molecule_kind_type, onfo_type, opbend_type, set_molecule_kind, & + shell_type, torsion_type, ub_type USE molecule_types_new, ONLY: get_molecule,& molecule_type USE pair_potential, ONLY: get_nonbond_storage,& @@ -89,12 +90,14 @@ MODULE force_fields_all force_field_unique_ub,& force_field_unique_tors,& force_field_unique_impr,& + force_field_unique_opbend,& force_field_unique_onfo,& force_field_pack_bond,& force_field_pack_bend,& force_field_pack_ub,& force_field_pack_tors,& force_field_pack_impr,& + force_field_pack_opbend,& force_field_pack_ch,& force_field_pack_shell,& force_field_pack_sbe,& @@ -645,10 +648,10 @@ CONTAINS ((name_atm_b)==(name_atm_b2)) .AND. & ((name_atm_c)==(name_atm_c2)) .AND. & ((name_atm_d)==(name_atm_d2))) .OR. & - (((name_atm_a)==(name_atm_d2)) .AND. & - ((name_atm_b)==(name_atm_c2)) .AND. & - ((name_atm_c)==(name_atm_b2)) .AND. & - ((name_atm_d)==(name_atm_a2))) ) THEN + (((name_atm_a)==(name_atm_a2)) .AND. & + ((name_atm_b)==(name_atm_b2)) .AND. & + ((name_atm_c)==(name_atm_d2)) .AND. & + ((name_atm_d)==(name_atm_c2))) ) THEN found = .TRUE. map_impr_kind(j) = map_impr_kind(k) EXIT @@ -675,6 +678,133 @@ CONTAINS END SUBROUTINE force_field_unique_impr +! ***************************************************************************** +!> \brief Determine the number of unique opbend kind and allocate opbend_kind_set +!> based on the present impropers. With each improper, there also +!> corresponds a opbend +! ***************************************************************************** + SUBROUTINE force_field_unique_opbend (particle_set, & + molecule_kind_set, molecule_set, ff_type, error) + + TYPE(particle_type), DIMENSION(:), & + POINTER :: particle_set + TYPE(molecule_kind_type), DIMENSION(:), & + POINTER :: molecule_kind_set + TYPE(molecule_type), DIMENSION(:), & + POINTER :: molecule_set + TYPE(force_field_type), INTENT(INOUT) :: ff_type + TYPE(cp_error_type), INTENT(INOUT) :: error + + CHARACTER(len=*), PARAMETER :: routineN = 'force_field_unique_opbend', & + routineP = moduleN//':'//routineN + + CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_a2, & + name_atm_b, name_atm_b2, name_atm_c, name_atm_c2, name_atm_d, & + name_atm_d2 + INTEGER :: atm_a, atm_b, atm_c, atm_d, & + counter, first, handle2, i, & + j, k, last, natom, nopbend, & + stat + INTEGER, DIMENSION(:), POINTER :: molecule_list + INTEGER, POINTER :: map_opbend_kind(:) + LOGICAL :: failure, found + TYPE(atomic_kind_type), POINTER :: atomic_kind + TYPE(molecule_kind_type), POINTER :: molecule_kind + TYPE(molecule_type), POINTER :: molecule + TYPE(opbend_kind_type), DIMENSION(:), & + POINTER :: opbend_kind_set + TYPE(opbend_type), DIMENSION(:), POINTER :: opbend_list + + failure = .FALSE. + + CALL timeset(routineN,handle2) + DO i=1,SIZE(molecule_kind_set) + molecule_kind => molecule_kind_set(i) + CALL get_molecule_kind(molecule_kind=molecule_kind,& + molecule_list=molecule_list,& + natom=natom,& + nopbend=nopbend,opbend_list=opbend_list) + molecule=>molecule_set(molecule_list(1)) + CALL get_molecule(molecule=molecule,first_atom=first,last_atom=last) + IF(nopbend>0) THEN + ALLOCATE(map_opbend_kind(nopbend),STAT=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + counter=0 + IF((ff_type%ff_type==do_ff_g96).OR.(ff_type%ff_type==do_ff_g87)) THEN + DO j=1,nopbend + map_opbend_kind(j)=j + END DO + counter=nopbend + ELSE + DO j=1,nopbend + atm_a = opbend_list(j)%a + atomic_kind => particle_set(atm_a+first-1)%atomic_kind + CALL get_atomic_kind(atomic_kind=atomic_kind,& + name=name_atm_a) + atm_b = opbend_list(j)%b + atomic_kind => particle_set(atm_b+first-1)%atomic_kind + CALL get_atomic_kind(atomic_kind=atomic_kind,& + name=name_atm_b) + atm_c = opbend_list(j)%c + atomic_kind => particle_set(atm_c+first-1)%atomic_kind + CALL get_atomic_kind(atomic_kind=atomic_kind,& + name=name_atm_c) + atm_d = opbend_list(j)%d + atomic_kind => particle_set(atm_d+first-1)%atomic_kind + CALL get_atomic_kind(atomic_kind=atomic_kind,& + name=name_atm_d) + found = .FALSE. + DO k=1,j-1 + atm_a = opbend_list(k)%a + atomic_kind => particle_set(atm_a+first-1)%atomic_kind + CALL get_atomic_kind(atomic_kind=atomic_kind,& + name=name_atm_a2) + atm_b = opbend_list(k)%b + atomic_kind => particle_set(atm_b+first-1)%atomic_kind + CALL get_atomic_kind(atomic_kind=atomic_kind,& + name=name_atm_b2) + atm_c = opbend_list(k)%c + atomic_kind => particle_set(atm_c+first-1)%atomic_kind + CALL get_atomic_kind(atomic_kind=atomic_kind,& + name=name_atm_c2) + atm_d = opbend_list(k)%d + atomic_kind => particle_set(atm_d+first-1)%atomic_kind + CALL get_atomic_kind(atomic_kind=atomic_kind,& + name=name_atm_d2) + IF( (((name_atm_a)==(name_atm_a2)) .AND. & + ((name_atm_b)==(name_atm_b2)) .AND. & + ((name_atm_c)==(name_atm_c2)) .AND. & + ((name_atm_d)==(name_atm_d2))) .OR. & + (((name_atm_a)==(name_atm_a2)) .AND. & + ((name_atm_b)==(name_atm_c2)) .AND. & + ((name_atm_c)==(name_atm_b2)) .AND. & + ((name_atm_d)==(name_atm_d2))) ) THEN + found = .TRUE. + map_opbend_kind(j) = map_opbend_kind(k) + EXIT + END IF + END DO + IF(.NOT.found) THEN + counter=counter+1 + map_opbend_kind(j) = counter + END IF + END DO + END IF + NULLIFY(opbend_kind_set) + CALL allocate_opbend_kind_set(opbend_kind_set,counter,error) + DO j=1,nopbend + opbend_list(j)%opbend_kind => opbend_kind_set(map_opbend_kind(j)) + END DO + CALL set_molecule_kind(molecule_kind=molecule_kind,& + opbend_kind_set=opbend_kind_set,opbend_list=opbend_list) + DEALLOCATE(map_opbend_kind,STAT=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + END IF + END DO + CALL timestop(handle2) + + END SUBROUTINE force_field_unique_opbend + ! ***************************************************************************** !> \brief Determine the number of unique onfo kind and allocate onfo_kind_set ! ***************************************************************************** @@ -1631,22 +1761,22 @@ CONTAINS IF(ASSOCIATED(inp_info%impr_a)) THEN DO k=1,SIZE(inp_info%impr_a) IF( ((inp_info%impr_a(k))==(name_atm_a)) .AND. & - ((inp_info%impr_d(k))==(name_atm_d)) .AND. & + ((inp_info%impr_b(k))==(name_atm_b)) .AND. & ((((inp_info%impr_c(k))==(name_atm_c)) .AND. & - ((inp_info%impr_b(k))==(name_atm_b))) .OR. & - (((inp_info%impr_c(k))==(name_atm_b)) .AND. & - ((inp_info%impr_b(k))==(name_atm_c)))) ) THEN + ((inp_info%impr_d(k))==(name_atm_d))) .OR. & + (((inp_info%impr_c(k))==(name_atm_d)) .AND. & + ((inp_info%impr_d(k))==(name_atm_c)))) ) THEN impr_list(j)%impr_kind%id_type = inp_info%impr_kind(k) impr_list(j)%impr_kind%k = inp_info%impr_k(k) IF ( ((inp_info%impr_c(k))==(name_atm_c)) .AND. & - ((inp_info%impr_b(k))==(name_atm_b)) ) THEN + ((inp_info%impr_d(k))==(name_atm_d)) ) THEN impr_list(j)%impr_kind%phi0 = inp_info%impr_phi0(k) ELSE impr_list(j)%impr_kind%phi0 = -inp_info%impr_phi0(k) ! alternative solutions: - ! - swap impr_list(j)%c with impr_list(j)%b and - ! name_atom_c with name_atom_b and - ! atm_c with atm_b + ! - swap impr_list(j)%c with impr_list(j)%d and + ! name_atom_c with name_atom_d and + ! atm_c with atm_d ! - introduce impr_list(j)%impr_kind%sign. if one, the ! sign of phi is not changed in mol_force.f90. if minus ! one, the sign of phi is changed in mol_force.f90 @@ -1665,7 +1795,7 @@ CONTAINS CALL store_FF_missing_par(atm1=TRIM(name_atm_a),& atm2=TRIM(name_atm_b),& atm3=TRIM(name_atm_c),& - atm4=TRIM(name_atm_c),& + atm4=TRIM(name_atm_d),& type_name="Improper",& array=Ainfo,& error=error) @@ -1689,6 +1819,138 @@ CONTAINS END SUBROUTINE force_field_pack_impr +! ***************************************************************************** +!> \brief Pack in opbend information needed for the force_field. +!> No loop over params for charmm, amber and gromos since these force +!> fields have no opbends +!> \author Louis Vanduyfhuys +! ***************************************************************************** + SUBROUTINE force_field_pack_opbend (particle_set, molecule_kind_set, molecule_set, & + Ainfo, inp_info, error) + + TYPE(particle_type), DIMENSION(:), & + POINTER :: particle_set + TYPE(molecule_kind_type), DIMENSION(:), & + POINTER :: molecule_kind_set + TYPE(molecule_type), DIMENSION(:), & + POINTER :: molecule_set + CHARACTER(LEN=default_string_length), & + DIMENSION(:), POINTER :: Ainfo + TYPE(input_info_type), POINTER :: inp_info + TYPE(cp_error_type), INTENT(INOUT) :: error + + CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_opbend', & + routineP = moduleN//':'//routineN + + CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_b, & + name_atm_c, name_atm_d + INTEGER :: atm_a, atm_b, atm_c, atm_d, & + first, handle2, i, j, k, & + last, natom, nopbend + INTEGER, DIMENSION(:), POINTER :: molecule_list + LOGICAL :: failure, found, only_qm + TYPE(atomic_kind_type), POINTER :: atomic_kind + TYPE(molecule_kind_type), POINTER :: molecule_kind + TYPE(molecule_type), POINTER :: molecule + TYPE(opbend_type), DIMENSION(:), POINTER :: opbend_list + + CALL timeset(routineN,handle2) + failure = .FALSE. + + DO i=1,SIZE(molecule_kind_set) + molecule_kind => molecule_kind_set(i) + CALL get_molecule_kind(molecule_kind=molecule_kind,& + molecule_list=molecule_list,& + natom=natom,& + nopbend=nopbend,opbend_list=opbend_list) + molecule=>molecule_set(molecule_list(1)) + + CALL get_molecule(molecule=molecule,first_atom=first,last_atom=last) + DO j=1,nopbend + atm_a = opbend_list(j)%a + atomic_kind => particle_set(atm_a+first-1)%atomic_kind + CALL get_atomic_kind(atomic_kind=atomic_kind,& + name=name_atm_a) + atm_b = opbend_list(j)%b + atomic_kind => particle_set(atm_b+first-1)%atomic_kind + CALL get_atomic_kind(atomic_kind=atomic_kind,& + name=name_atm_b) + atm_c = opbend_list(j)%c + atomic_kind => particle_set(atm_c+first-1)%atomic_kind + CALL get_atomic_kind(atomic_kind=atomic_kind,& + name=name_atm_c) + atm_d = opbend_list(j)%d + atomic_kind => particle_set(atm_d+first-1)%atomic_kind + CALL get_atomic_kind(atomic_kind=atomic_kind,& + name=name_atm_d) + found = .FALSE. + only_qm = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b, id3=name_atm_c, id4=name_atm_d) + CALL uppercase(name_atm_a) + CALL uppercase(name_atm_b) + CALL uppercase(name_atm_c) + CALL uppercase(name_atm_d) + + ! always have the input param last to overwrite all the other ones + IF(ASSOCIATED(inp_info%opbend_a)) THEN + DO k=1,SIZE(inp_info%opbend_a) + IF( ((inp_info%opbend_a(k))==(name_atm_a)) .AND. & + ((inp_info%opbend_d(k))==(name_atm_d)) .AND. & + ((((inp_info%opbend_c(k))==(name_atm_c)) .AND. & + ((inp_info%opbend_b(k))==(name_atm_b))) .OR. & + (((inp_info%opbend_c(k))==(name_atm_b)) .AND. & + ((inp_info%opbend_b(k))==(name_atm_c)))) ) THEN + opbend_list(j)%opbend_kind%id_type = inp_info%opbend_kind(k) + opbend_list(j)%opbend_kind%k = inp_info%opbend_k(k) + IF ( ((inp_info%opbend_c(k))==(name_atm_c)) .AND. & + ((inp_info%opbend_b(k))==(name_atm_b)) ) THEN + opbend_list(j)%opbend_kind%phi0 = inp_info%opbend_phi0(k) + ELSE + opbend_list(j)%opbend_kind%phi0 = -inp_info%opbend_phi0(k) + ! alternative solutions: + ! - swap opbend_list(j)%c with opbend_list(j)%b and + ! name_atom_c with name_atom_b and + ! atm_c with atm_b + ! - introduce opbend_list(j)%opbend_kind%sign. if one, the + ! sign of phi is not changed in mol_force.f90. if minus + ! one, the sign of phi is changed in mol_force.f90 + + END IF + CALL issue_duplications(found,routineP, "Out of plane bend", name_atm_a, name_atm_b,& + name_atm_c,name_atm_d, error=error) + found = .TRUE. + EXIT + END IF + END DO + END IF + + IF(.NOT.found) THEN + CALL store_FF_missing_par(atm1=TRIM(name_atm_a),& + atm2=TRIM(name_atm_b),& + atm3=TRIM(name_atm_c),& + atm4=TRIM(name_atm_d),& + type_name="Out of plane bend",& + array=Ainfo,& + error=error) + opbend_list(j)%opbend_kind%k = 0.0_dp + opbend_list(j)%opbend_kind%phi0 = 0.0_dp + opbend_list(j)%opbend_kind%id_type = do_ff_undef + opbend_list(j)%id_type = do_ff_undef + END IF + ! + ! QM/MM modifications + ! + IF (only_qm) THEN + opbend_list(j)%opbend_kind%id_type = do_ff_undef + opbend_list(j)%id_type = do_ff_undef + END IF + + END DO + CALL set_molecule_kind(molecule_kind=molecule_kind,opbend_list=opbend_list) + END DO + CALL timestop(handle2) + + END SUBROUTINE force_field_pack_opbend + ! ***************************************************************************** !> \brief Set up atomic_kind_set()%elp_potentail%qeff (PART 2) and shell potential !> parameters @@ -3001,6 +3263,21 @@ CPSourceFileRef,& ((atm1==my_atm1).AND.(atm2==my_atm4).AND.(atm3==my_atm2).AND.(atm4==my_atm3)).OR.& ((atm1==my_atm1).AND.(atm2==my_atm2).AND.(atm3==my_atm4).AND.(atm4==my_atm3)))& found = .TRUE. + + CASE("Out of plane bend") + IF (INDEX(array(i)(21:39),"Out of plane bend")==0) CYCLE + my_atm1 = array(i)(41:44) + my_atm2 = array(i)(46:49) + my_atm3 = array(i)(51:54) + my_atm4 = array(i)(56:59) + CALL compress(my_atm1,.TRUE.) + CALL compress(my_atm2,.TRUE.) + CALL compress(my_atm3,.TRUE.) + CALL compress(my_atm4,.TRUE.) + IF (((atm1==my_atm1).AND.(atm2==my_atm2).AND.(atm3==my_atm3).AND.(atm4==my_atm4)).OR.& + ((atm1==my_atm1).AND.(atm2==my_atm3).AND.(atm3==my_atm2).AND.(atm4==my_atm4)))& + found = .TRUE. + CASE("Charge") IF (INDEX(array(i)(21:39),"Charge")==0) CYCLE my_atm1 = array(i)(41:44) diff --git a/src/force_fields_input.F b/src/force_fields_input.F index 2a54241271..8584ef911f 100644 --- a/src/force_fields_input.F +++ b/src/force_fields_input.F @@ -87,7 +87,7 @@ CONTAINS routineP = moduleN//':'//routineN INTEGER :: nb4, nbends, nbm, nbmhft, nbonds, nchg, neam, ngd, ngp, nimpr, & - nipbv, nlj, nshell, ntersoff, ntors, ntot, nubs, nwl + nipbv, nlj, nopbend, nshell, ntersoff, ntors, ntot, nubs, nwl LOGICAL :: explicit, failure, & unique_spline REAL(KIND=dp) :: min_eps_spline_allowed @@ -365,6 +365,22 @@ CPSourceFileRef,& inp_info%impr_c,inp_info%impr_d,inp_info%impr_k,inp_info%impr_phi0,& tmp_section, ntot, error ) END IF + + tmp_section => section_vals_get_subs_vals(ff_section,"OPBEND",error=error) + CALL section_vals_get(tmp_section,explicit=explicit,n_repetition=nopbend,error=error) + IF (explicit) THEN + ntot=0 + CALL reallocate(inp_info%opbend_kind,1,nopbend) + CALL reallocate(inp_info%opbend_a,1,nopbend) + CALL reallocate(inp_info%opbend_b,1,nopbend) + CALL reallocate(inp_info%opbend_c,1,nopbend) + CALL reallocate(inp_info%opbend_d,1,nopbend) + CALL reallocate(inp_info%opbend_k,1,nopbend) + CALL reallocate(inp_info%opbend_phi0,1,nopbend) + CALL read_opbend_section(inp_info%opbend_kind,inp_info%opbend_a, inp_info%opbend_b,& + inp_info%opbend_c,inp_info%opbend_d,inp_info%opbend_k,inp_info%opbend_phi0,& + tmp_section, ntot, error ) + END IF END IF END SUBROUTINE read_force_field_section1 @@ -1341,6 +1357,48 @@ CPSourceFileRef,& END DO END SUBROUTINE read_improper_section +! ***************************************************************************** +!> \brief Reads the OPBEND section +!> \author louis vanduyfhuys +! ***************************************************************************** + SUBROUTINE read_opbend_section(opbend_kind,opbend_a, opbend_b, opbend_c, opbend_d, opbend_k,& + opbend_phi0, section, start, error ) + INTEGER, DIMENSION(:), POINTER :: opbend_kind + CHARACTER(LEN=default_string_length), & + DIMENSION(:), POINTER :: opbend_a, opbend_b, opbend_c, & + opbend_d + REAL(KIND=dp), DIMENSION(:), POINTER :: opbend_k, opbend_phi0 + TYPE(section_vals_type), POINTER :: section + INTEGER, INTENT(IN) :: start + TYPE(cp_error_type), INTENT(inout) :: error + + CHARACTER(len=*), PARAMETER :: routineN = 'read_opbend_section', & + routineP = moduleN//':'//routineN + + CHARACTER(LEN=default_string_length), & + DIMENSION(:), POINTER :: atm_names + INTEGER :: isec, n_items + LOGICAL :: failure + + failure = .FALSE. + NULLIFY( atm_names) + CALL section_vals_get(section,n_repetition=n_items,error=error) + DO isec = 1, n_items + CALL section_vals_val_get(section,"KIND",i_rep_section=isec,i_val=opbend_kind(start+isec),error=error) + CALL section_vals_val_get(section,"ATOMS",i_rep_section=isec,c_vals=atm_names,error=error) + opbend_a(start+isec) = atm_names(1) + opbend_b(start+isec) = atm_names(2) + opbend_c(start+isec) = atm_names(3) + opbend_d(start+isec) = atm_names(4) + CALL uppercase(opbend_a(start+isec)) + CALL uppercase(opbend_b(start+isec)) + CALL uppercase(opbend_c(start+isec)) + CALL uppercase(opbend_d(start+isec)) + CALL section_vals_val_get(section,"K",i_rep_section=isec,r_val=opbend_k(start+isec),error=error) + CALL section_vals_val_get(section,"PHI0",i_rep_section=isec,r_val=opbend_phi0(start+isec),error=error) + END DO + END SUBROUTINE read_opbend_section + ! ***************************************************************************** !> \brief Reads the force_field input section !> \par History diff --git a/src/force_fields_util.F b/src/force_fields_util.F index fcfe945df2..a74377fc72 100644 --- a/src/force_fields_util.F +++ b/src/force_fields_util.F @@ -29,20 +29,22 @@ MODULE force_fields_util fist_nonbond_env_type USE force_field_types, ONLY: & allocate_bend_kind_set, allocate_bond_kind_set, & - allocate_impr_kind_set, allocate_torsion_kind_set, & - allocate_ub_kind_set, amber_info_type, bend_kind_type, bond_kind_type, & - charmm_info_type, deallocate_bend_kind_set, deallocate_bond_kind_set, & - force_field_type, gromos_info_type, impr_kind_dealloc_ref, & - impr_kind_type, input_info_type, torsion_kind_dealloc_ref, & + allocate_impr_kind_set, allocate_opbend_kind_set, & + allocate_torsion_kind_set, allocate_ub_kind_set, amber_info_type, & + bend_kind_type, bond_kind_type, charmm_info_type, & + deallocate_bend_kind_set, deallocate_bond_kind_set, force_field_type, & + gromos_info_type, impr_kind_dealloc_ref, impr_kind_type, & + input_info_type, opbend_kind_type, torsion_kind_dealloc_ref, & torsion_kind_type, ub_kind_dealloc_ref, ub_kind_type USE force_fields_all, ONLY: & force_field_pack_bend, force_field_pack_bond, force_field_pack_ch, & - force_field_pack_impr, force_field_pack_onfo, force_field_pack_pp, & - force_field_pack_ppse, force_field_pack_sbe, force_field_pack_shell, & - force_field_pack_sne, force_field_pack_tors, force_field_pack_ub, & - force_field_unique_bend, force_field_unique_bond, & + force_field_pack_impr, force_field_pack_onfo, force_field_pack_opbend, & + force_field_pack_pp, force_field_pack_ppse, force_field_pack_sbe, & + force_field_pack_shell, force_field_pack_sne, force_field_pack_tors, & + force_field_pack_ub, force_field_unique_bend, force_field_unique_bond, & force_field_unique_impr, force_field_unique_onfo, & - force_field_unique_tors, force_field_unique_ub + force_field_unique_opbend, force_field_unique_tors, & + force_field_unique_ub USE input_constants, ONLY: do_ff_undef USE input_section_types, ONLY: section_vals_type,& section_vals_val_get @@ -53,7 +55,8 @@ MODULE force_fields_util USE molecule_kind_types, ONLY: & atom_type, bend_type, bond_type, colvar_constraint_type, & g3x3_constraint_type, g4x6_constraint_type, get_molecule_kind, & - impr_type, molecule_kind_type, set_molecule_kind, torsion_type, ub_type + impr_type, molecule_kind_type, opbend_type, set_molecule_kind, & + torsion_type, ub_type USE molecule_types_new, ONLY: molecule_type USE pair_potential_types, ONLY: pair_potential_pp_type USE particle_types, ONLY: particle_type @@ -185,48 +188,61 @@ CONTAINS ff_type, error) !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- - ! 6. Determine the number of unique onfo kind and allocate onfo_kind_set + ! 6. Determine the number of unique opbend kind and allocate opbend_kind_set + !----------------------------------------------------------------------------- + CALL force_field_unique_opbend (particle_set, molecule_kind_set, molecule_set, & + ff_type, error) + !----------------------------------------------------------------------------- + !----------------------------------------------------------------------------- + ! 7. Determine the number of unique onfo kind and allocate onfo_kind_set !----------------------------------------------------------------------------- CALL force_field_unique_onfo (particle_set, molecule_kind_set,& molecule_set, error) !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- - ! 7. Bonds + ! 8. Bonds !----------------------------------------------------------------------------- CALL force_field_pack_bond (particle_set, molecule_kind_set, molecule_set, & fatal, Ainfo, chm_info, inp_info, gro_info, amb_info, error) !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- - ! 8. Bends + ! 9. Bends !----------------------------------------------------------------------------- CALL force_field_pack_bend (particle_set, molecule_kind_set, molecule_set, & fatal, Ainfo, chm_info, inp_info, gro_info, amb_info, error) ! Give information and abort if any bond or angle parameter is missing.. - CALL release_FF_missing_par(fatal,AInfo,output_unit,iw,error) + CALL release_FF_missing_par(fatal,AInfo,output_unit,iw,error) !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- - ! 9. Urey-Bradley + ! 10. Urey-Bradley !----------------------------------------------------------------------------- CALL force_field_pack_ub (particle_set, molecule_kind_set, molecule_set, & Ainfo, chm_info, inp_info, iw, error) !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- - ! 10. Torsions + ! 11. Torsions !----------------------------------------------------------------------------- CALL force_field_pack_tors (particle_set, molecule_kind_set, molecule_set, & Ainfo, chm_info, inp_info, gro_info, amb_info, iw, error) !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- - ! 11. Impropers + ! 12. Impropers !----------------------------------------------------------------------------- CALL force_field_pack_impr (particle_set, molecule_kind_set, molecule_set, & Ainfo, chm_info, inp_info, gro_info, amb_info, error) - ! Give information only if any Urey-Bradley, Torsion or improper is missing + !----------------------------------------------------------------------------- + !----------------------------------------------------------------------------- + ! 13. Out of plane bends + !----------------------------------------------------------------------------- + CALL force_field_pack_opbend (particle_set, molecule_kind_set, molecule_set, & + Ainfo, inp_info, error) + + ! Give information only if any Urey-Bradley, Torsion, improper or opbend is missing ! continue calculation.. CALL release_FF_missing_par(fatal,AInfo,output_unit,iw,error) !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- - ! 12. Set up atomic_kind_set()%elp_potentail%qeff (PART 2) + ! 14. Set up atomic_kind_set()%elp_potentail%qeff (PART 2) !----------------------------------------------------------------------------- CALL force_field_pack_ch (atomic_kind_set, qmmm_env, fatal, iw, iw4,& Ainfo, my_qmmm, inp_info, error) @@ -235,14 +251,14 @@ CONTAINS CALL release_FF_missing_par(fatal,AInfo,output_unit,iw,error) !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- - ! 13. Set up Shell potential + ! 15. Set up Shell potential !----------------------------------------------------------------------------- CALL force_field_pack_shell ( particle_set, atomic_kind_set,& molecule_kind_set, molecule_set, root_section, subsys_section,& shell_particle_set, core_particle_set, cell, iw, inp_info, error) !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- - ! 14. Set spline_bond_env, use for ONFO interaction + ! 16. Set spline_bond_env, use for ONFO interaction !----------------------------------------------------------------------------- IF (ff_type%nonbonded14) THEN NULLIFY(spline_bond_env) @@ -253,14 +269,14 @@ CONTAINS CALL release_FF_missing_par(fatal,AInfo,output_unit,iw,error) !-------------------------------------------------------------------------- !-------------------------------------------------------------------------- - ! 14.(a) ONFO + ! 16.(a) ONFO !-------------------------------------------------------------------------- CALL force_field_pack_onfo(particle_set, molecule_kind_set, molecule_set,& iw, potparm_bond, error) END IF !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- - ! 15. Set spline_nonbond_env + ! 17. Set spline_nonbond_env !----------------------------------------------------------------------------- IF (ff_type%nonbonded) THEN CALL force_field_pack_sne (atomic_kind_set, ff_type, qmmm_env, ewald_env, & @@ -269,7 +285,7 @@ CONTAINS ! Give information and abort if any pair potential spline is missing.. CALL release_FF_missing_par(fatal,AInfo,output_unit,iw,error) !-------------------------------------------------------------------------- - ! 15.(a) create the pair potential spline environment + ! 17.(a) create the pair potential spline environment !-------------------------------------------------------------------------- NULLIFY(spline_nonbond_env, spline_c_coul_env, spline_cs_coul_env, & spline_sc_coul_env) @@ -278,7 +294,7 @@ CONTAINS iw4, spline_nonbond_env, spline_c_coul_env, spline_cs_coul_env,& spline_sc_coul_env, potparm_nonbond, do_zbl, error) !-------------------------------------------------------------------------- - ! 15.(b) Setup of the pair potential spline environment (Evtl. Core/Shell) + ! 17.(b) Setup of the pair potential spline environment (Evtl. Core/Shell) !-------------------------------------------------------------------------- CALL force_field_pack_pp (atomic_kind_set, spline_nonbond_env, spline_c_coul_env,& spline_cs_coul_env, spline_sc_coul_env, potparm_nonbond, shell_present,& @@ -334,7 +350,7 @@ CONTAINS WRITE(output_unit,*) WRITE(output_unit,'(T2,"FORCEFIELD|",A)')& " WARNING: A non Critical ForceField parameter is missing! CP2K GOES!",& - " Non critical parameters are:Urey-Bradley,Torsions,Impropers and 1-4",& + " Non critical parameters are:Urey-Bradley,Torsions,Impropers, Opbends and 1-4",& " All missing parameters will not contribute to the potential energy!" IF (fatal.OR.iw>0) THEN WRITE(output_unit,*) @@ -396,7 +412,7 @@ CONTAINS qeff_mol = 0.0_dp qeff_sum = 0.0_dp mass_sum = 0.0_dp - + !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- ! 1. Sum of [qeff,mass] for each molecule_kind @@ -418,7 +434,7 @@ CONTAINS CALL set_molecule_kind(molecule_kind=molecule_kind,charge=qeff_mol,mass=mass_mol) IF(iw>0) WRITE(iw,*) " Mol Kind ",TRIM(molname)," charge = ",qeff_mol END DO - + !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- ! 2. Sum of [qeff,mass] for particle_set @@ -434,7 +450,7 @@ CONTAINS END DO IF(iw>0) WRITE(iw,*) " Total system charge = ",qeff_sum IF(iw>0) WRITE(iw,*) " Total system mass = ",mass_sum - + CALL cp_print_key_finished_output(iw,logger,mm_section,& "PRINT%FF_INFO",error=error) CALL timestop(handle) @@ -454,9 +470,9 @@ CONTAINS routineP = moduleN//':'//routineN INTEGER :: atm2_a, atm2_b, atm2_c, atm_a, atm_b, atm_c, atm_d, counter, & - handle, i, ibend, ibond, icolv, ig3x3, ig4x6, iimpr, ikind, itorsion, & - iub, iw, j, k, nbend, nbond, newkind, ng3x3, ng4x6, nimpr, ntorsion, & - nub, stat + handle, i, ibend, ibond, icolv, ig3x3, ig4x6, iimpr, ikind, iopbend, & + itorsion, iub, iw, j, k, nbend, nbond, newkind, ng3x3, ng4x6, nimpr, & + nopbend, ntorsion, nub, stat INTEGER, POINTER :: bad1(:), bad2(:) LOGICAL :: failure, unsetme, valid_kind TYPE(bend_kind_type), DIMENSION(:), & @@ -479,6 +495,10 @@ CONTAINS new_impr_kind_set TYPE(impr_type), DIMENSION(:), POINTER :: impr_list, new_impr_list TYPE(molecule_kind_type), POINTER :: molecule_kind + TYPE(opbend_kind_type), DIMENSION(:), & + POINTER :: new_opbend_kind_set, & + opbend_kind_set + TYPE(opbend_type), DIMENSION(:), POINTER :: new_opbend_list, opbend_list TYPE(torsion_kind_type), DIMENSION(:), & POINTER :: new_torsion_kind_set, & torsion_kind_set @@ -1130,6 +1150,93 @@ CONTAINS END IF END DO + !----------------------------------------------------------------------------- + !----------------------------------------------------------------------------- + ! 11. Count the number of UNSET opbends kinds there are + !----------------------------------------------------------------------------- + DO ikind=1,SIZE(molecule_kind_set) + molecule_kind => molecule_kind_set(ikind) + CALL get_molecule_kind(molecule_kind=molecule_kind,& + nopbend=nopbend,& + opbend_kind_set=opbend_kind_set,& + opbend_list=opbend_list) + IF(nopbend>0) THEN + IF(iw>0) WRITE(iw,*) " Mol(",ikind,") Old OPBEND Count: ",& + SIZE(opbend_list),SIZE(opbend_kind_set) + NULLIFY(bad1,bad2) + ALLOCATE(bad1(SIZE(opbend_kind_set)),STAT=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + bad1(:)=0 + DO iopbend=1,SIZE(opbend_kind_set) + unsetme=.FALSE. + IF(opbend_kind_set(iopbend)%id_type==do_ff_undef) unsetme = .TRUE. + valid_kind=.FALSE. + DO i=1,SIZE(opbend_list) + IF(opbend_list(i)%id_type/=do_ff_undef.AND.& + opbend_list(i)%opbend_kind%kind_number==iopbend) THEN + valid_kind=.TRUE. + EXIT + END IF + END DO + IF(.NOT.valid_kind) unsetme = .TRUE. + IF(unsetme) bad1(iopbend) = 1 + END DO + IF(SUM(bad1)/=0) THEN + counter = SIZE(opbend_kind_set)-SUM(bad1) + CALL allocate_opbend_kind_set(new_opbend_kind_set,counter,error) + counter=0 + DO iopbend=1,SIZE(opbend_kind_set) + IF(bad1(iopbend)==0) THEN + counter=counter+1 + new_opbend_kind_set(counter)= opbend_kind_set(iopbend) + END IF + END DO + counter=0 + ALLOCATE(bad2(SIZE(opbend_list)),STAT=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + bad2(:)=0 + DO iopbend=1,SIZE(opbend_list) + unsetme = .FALSE. + IF(opbend_list(iopbend)%opbend_kind%id_type==do_ff_undef) unsetme = .TRUE. + IF(opbend_list(iopbend)%id_type==do_ff_undef) unsetme = .TRUE. + IF(unsetme) bad2(iopbend) = 1 + END DO + IF(SUM(bad2)/=0) THEN + counter = SIZE(opbend_list)-SUM(bad2) + ALLOCATE(new_opbend_list(counter),STAT=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + counter=0 + DO iopbend=1,SIZE(opbend_list) + IF(bad2(iopbend)==0) THEN + counter=counter+1 + new_opbend_list(counter) = opbend_list(iopbend) + newkind = opbend_list(iopbend)%opbend_kind%kind_number + newkind = newkind - SUM(bad1(1:newkind)) + new_opbend_list(counter)%opbend_kind => new_opbend_kind_set(newkind) + END IF + END DO + CALL set_molecule_kind(molecule_kind=molecule_kind,& + nopbend=SIZE(new_opbend_list),& + opbend_kind_set=new_opbend_kind_set,& + opbend_list=new_opbend_list) + DO iopbend=1,SIZE(new_opbend_kind_set) + new_opbend_kind_set(iopbend)%kind_number=iopbend + END DO + END IF + DEALLOCATE(bad2,STAT=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + DEALLOCATE(opbend_kind_set,STAT=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + DEALLOCATE(opbend_list,STAT=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + IF(iw>0) WRITE(iw,*) " Mol(",ikind,") New OPBEND Count: ",& + SIZE(new_opbend_list),SIZE(new_opbend_kind_set) + END IF + DEALLOCATE(bad1,STAT=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + END IF + END DO + !--------------------------------------------------------------------------- !--------------------------------------------------------------------------- ! 12. Count the number of UNSET NONBOND14 kinds there are @@ -1191,7 +1298,7 @@ CONTAINS ELSE CALL section_vals_val_get(gen_section,"VARIABLES",i_rep_section=irep,c_vals=my_var,error=error) END IF - IF (ASSOCIATED(my_var)) THEN + IF (ASSOCIATED(my_var)) THEN mydim = SIZE(my_var) END IF IF (PRESENT(size_variables)) THEN @@ -1243,7 +1350,7 @@ CONTAINS IF (mydim>0) THEN parameters(1:SIZE(my_var)) = my_var values(1:SIZE(my_var)) = 0.0_dp - IF (PRESENT(var_values)) THEN + IF (PRESENT(var_values)) THEN CPPrecondition(SIZE(var_values)==SIZE(my_var),cp_failure_level,routineP,error,failure) values(1:SIZE(my_var)) = var_values END IF diff --git a/src/input_constants.F b/src/input_constants.F index 322366bbaa..9618baa34c 100644 --- a/src/input_constants.F +++ b/src/input_constants.F @@ -8,7 +8,7 @@ !> circular dependencies !> \par History !> 02.2005 created [fawzi&Joost] -!> 10.2008 Teodoro Laino [tlaino] - University of Zurich +!> 10.2008 Teodoro Laino [tlaino] - University of Zurich !> Cleaned this file of all arrays of strings. Here must be kept !> only integer/real global constants ! ***************************************************************************** @@ -44,7 +44,7 @@ MODULE input_constants do_ep=6,& do_mixed=7 - ! QMMM + ! QMMM REAL(KIND=dp), PARAMETER, PUBLIC :: RADIUS_QMMM_DEFAULT = 0.80_dp,& ! Angstrom ALPHA_IMOMM_DEFAULT = 1.38_dp,& CHARGE_SCALE_FACTOR = 1.00_dp @@ -58,31 +58,31 @@ MODULE input_constants INTEGER, PARAMETER, PUBLIC :: weight_type_unit=0,& weight_type_mass=1 - ! Blacs type of distribution + ! Blacs type of distribution INTEGER, PARAMETER, PUBLIC :: BLACS_GRID_SQUARE = 1,& BLACS_GRID_ROW = 2,& BLACS_GRID_COL = 3 - ! EWALD methods + ! EWALD methods INTEGER, PARAMETER, PUBLIC :: do_ewald_none=1,& do_ewald_ewald=2,& do_ewald_pme=3,& do_ewald_spme=4 - + INTEGER, PARAMETER, PUBLIC :: do_constr_none=1,& do_constr_atomic=2,& do_constr_molec=3 - + INTEGER, PARAMETER, PUBLIC :: do_cell_cp2k=1,& do_cell_xsc=2 INTEGER, PARAMETER, PUBLIC :: plane_def_atoms=0,& - plane_def_vec=1 + plane_def_vec=1 INTEGER, PARAMETER, PUBLIC :: do_rep_interleaved=0,& - do_rep_blocked=1 + do_rep_blocked=1 - ! MD Ensemble + ! MD Ensemble INTEGER, PARAMETER, PUBLIC :: nve_ensemble=1,& nvt_ensemble=2,& npt_i_ensemble=3,& @@ -95,14 +95,14 @@ MODULE input_constants langevin_ensemble=10,& npe_f_ensemble=11,& npe_i_ensemble=12 - + INTEGER, PARAMETER, PUBLIC :: do_thermo_same_as_part=100,& do_thermo_nose=101,& do_thermo_csvr=102,& do_thermo_gle=103 - - ! Run Types + + ! Run Types INTEGER, PARAMETER, PUBLIC :: none_run=1,& energy_run=2,& energy_force_run=3,& @@ -119,8 +119,8 @@ MODULE input_constants cell_opt_run=14,& real_time_propagation=15,& ehrenfest=16 - - ! Run Types of Atom Code + + ! Run Types of Atom Code INTEGER, PARAMETER, PUBLIC :: atom_no_run=1,& atom_energy_run=2,& atom_basis_run=3,& @@ -238,10 +238,10 @@ MODULE input_constants INTEGER, PARAMETER, PUBLIC :: rsgrid_distributed=0,& rsgrid_replicated=1,& rsgrid_automatic=2 - + INTEGER, PARAMETER, PUBLIC :: do_par_atom=0,& do_par_grid=1 - ! Thermostats + ! Thermostats INTEGER, PARAMETER, PUBLIC :: do_thermo_no_communication=0,& do_thermo_communication=1,& do_thermo_only_master=2 @@ -249,25 +249,25 @@ MODULE input_constants ! pint coordinate transformations INTEGER, PARAMETER, PUBLIC :: transformation_normal=1,& transformation_stage=2 - - ! Free Energy methods + + ! Free Energy methods INTEGER, PARAMETER, PUBLIC :: do_fe_meta=0,& do_fe_ui=1,& do_fe_ac=2 - + INTEGER, PARAMETER, PUBLIC :: do_wall_none=0,& do_wall_reflective=1,& do_wall_quadratic=2,& do_wall_quartic=3,& do_wall_gaussian=4 - + INTEGER, PARAMETER, PUBLIC :: do_wall_p=0,& do_wall_m=1 - ! Output formats + ! Output formats INTEGER, PARAMETER, PUBLIC :: dump_atomic=1,& dump_xmol=2,& dump_dcd=3 - + INTEGER, PARAMETER, PUBLIC :: do_qmmm_none=0,& do_qmmm_coulomb=1,& do_qmmm_gauss=2,& @@ -280,7 +280,7 @@ MODULE input_constants do_qmmm_wall_none=0,& do_qmmm_wall_reflective=1,& do_qmmm_wall_quadratic=2 - + INTEGER, PARAMETER, PUBLIC :: do_region_global=0, & do_region_molecule=1,& do_region_massive=2,& @@ -305,7 +305,7 @@ MODULE input_constants do_coord_crd=4,& do_coord_cif=5,& do_coord_xtl=6 - + INTEGER, PARAMETER, PUBLIC :: do_ff_undef=0,& do_ff_quartic=1,& do_ff_g96=2,& @@ -316,30 +316,32 @@ MODULE input_constants do_ff_cubic=7,& do_ff_mixed_bend_stretch=8,& do_ff_amber=9,& - do_ff_mm3=10 - + do_ff_mm2=10,& + do_ff_mm3=11,& + do_ff_mm4=12 + INTEGER, PARAMETER, PUBLIC :: do_stress_none=0,& do_stress_analytical=1,& do_stress_numerical=2,& do_stress_diagonal_anal=3,& do_stress_diagonal_numer=4 - + INTEGER, PARAMETER, PUBLIC :: do_mc_traditional=0,& do_mc_gemc_nvt=1,& do_mc_gemc_npt=2,& do_mc_virial=3 - + INTEGER, PARAMETER, PUBLIC :: do_multipole_none=-1,& do_multipole_charge=0,& do_multipole_dipole=1,& do_multipole_quadrupole=2 - + INTEGER, PARAMETER, PUBLIC :: do_loc_none=0,& do_loc_jacobi=1,& do_loc_crazy=2,& do_loc_direct=3,& do_loc_l1_norm_sd=4 - + INTEGER, PARAMETER, PUBLIC :: do_loc_min=0,& do_loc_max=1,& op_loc_berry=0,& @@ -361,12 +363,12 @@ MODULE input_constants orb_dx2=7,& orb_dy2=8,& orb_dz2=9 - + INTEGER, PARAMETER, PUBLIC :: use_mom_ref_com=0,& use_mom_ref_coac=1,& use_mom_ref_user=2,& use_mom_ref_zero=3 - + INTEGER, PARAMETER, PUBLIC :: current_gauge_r=1,& current_gauge_r_and_step_func=2,& current_gauge_atom=3,& @@ -374,7 +376,7 @@ MODULE input_constants current_orb_center_common=2,& current_orb_center_atom=3,& current_orb_center_box=4 - + INTEGER, PARAMETER, PUBLIC :: ot_mini_sd=1,ot_mini_cg=2,ot_mini_diis=3 INTEGER, PARAMETER, PUBLIC :: ot_algo_taylor_or_diag=1,ot_algo_irac=2 INTEGER, PARAMETER, PUBLIC :: ot_chol_irac=1,ot_poly_irac=2,ot_lwdn_irac=3 @@ -392,7 +394,7 @@ MODULE input_constants ot_precond_solver_inv_chol=2,& ot_precond_solver_direct=3,& ot_precond_solver_sainv=4 - + INTEGER, PARAMETER, PUBLIC :: diag_standard=1,& diag_ot=2, diag_block_krylov=3 @@ -404,22 +406,22 @@ MODULE input_constants INTEGER, PUBLIC, PARAMETER :: default_lbfgs_method_id=1,& default_bfgs_method_id=2,& default_cg_method_id=3 - + INTEGER, PUBLIC, PARAMETER :: default_minimization_method_id=1,& default_ts_method_id=2,& default_cell_method_id=3,& default_shellcore_method_id=4 - + INTEGER, PUBLIC, PARAMETER :: none_ts_method_id=0,& default_dimer_method_id=1 - + INTEGER, PUBLIC, PARAMETER :: do_first_rotation_step=0,& do_second_rotation_step=1,& do_third_rotation_step=2 - + INTEGER, PUBLIC, PARAMETER :: default_cell_geo_opt_id=0,& default_cell_md_id=1 - + INTEGER, PARAMETER, PUBLIC :: add_last_no=0,& add_last_numeric=1,& add_last_symbolic=2 @@ -428,7 +430,7 @@ MODULE input_constants medium_print_level=2,& high_print_level=3,& debug_print_level=4 - + INTEGER, PARAMETER, PUBLIC :: xc_rho_no_smooth=0,& xc_rho_spline2_smooth=1,& xc_rho_spline3_smooth=2,& @@ -464,7 +466,7 @@ MODULE input_constants sic_ad=3,& sic_eo=4 INTEGER, PARAMETER, PUBLIC :: sic_list_all=1,& - sic_list_unpaired=2 + sic_list_unpaired=2 INTEGER, PARAMETER, PUBLIC :: tddfpt_singlet=0,& tddfpt_triplet=1 INTEGER, PARAMETER, PUBLIC :: tddfpt_lanczos=0,& @@ -479,7 +481,7 @@ MODULE input_constants INTEGER, PARAMETER, PUBLIC :: xc_vdw_fun_none=100,& xc_vdw_fun_pairpot=101 INTEGER, PARAMETER, PUBLIC :: vdw_pairpot_grimme=1 - + ! Input constants for relativistic calculations INTEGER, PARAMETER, PUBLIC :: rel_none=0,& rel_dkh=1 @@ -488,17 +490,17 @@ MODULE input_constants rel_trans_atom=2 INTEGER, PARAMETER, PUBLIC :: rel_pot_full=0,& rel_pot_erfc=1 - + INTEGER, PARAMETER, PUBLIC :: do_add=1,& do_remove=2 - + INTEGER, PARAMETER, PUBLIC :: no_precond=0,& precond_spl3_aint=1,& precond_spl3_1=2,& precond_spl3_aint2=3,& precond_spl3_2=4,& precond_spl3_3=5 - + INTEGER, PARAMETER, PUBLIC :: PERIODIC3D = 1000,& ANALYTIC2D = 1001,& ANALYTIC1D = 1002,& @@ -516,7 +518,7 @@ MODULE input_constants WAVELET2D = 1401,& WAVELET1D = 1402,& WAVELET0D = 1403 - + INTEGER, PARAMETER, PUBLIC :: use_none=0,& use_periodic = 1,& use_analytic = 2,& @@ -524,10 +526,10 @@ MODULE input_constants use_hockney=5,& use_multipole=4,& use_wavelet=6 - + INTEGER, PARAMETER, PUBLIC :: use_rs_grid=0,& use_gs_grid = 1 - + INTEGER, PARAMETER, PUBLIC :: xgga_b88=101,& xgga_b88x=102,& xgga_pw86=103,& @@ -556,7 +558,7 @@ MODULE input_constants c_pzvmc=214,& c_vwn=215,& cgga_pbe=216 - + INTEGER, PARAMETER, PUBLIC :: ke_ol1=301,& ke_ol2=302,& ke_llp=303,& @@ -568,18 +570,18 @@ MODULE input_constants ke_revpbe=309,& ke_tf=310,& ke_tfw=311 - + INTEGER, PARAMETER, PUBLIC :: pw_orig=209,& pw_dmc=210,& pw_vmc=211 - + INTEGER, PARAMETER, PUBLIC :: pz_orig=212,& pz_dmc=213,& pz_vmc=214 - + INTEGER, PARAMETER, PUBLIC :: xc_b97_orig=401,& xc_b97_grimme=402 - + INTEGER, PARAMETER, PUBLIC :: use_perd_x=0,& use_perd_y=1,& use_perd_z=2,& @@ -588,7 +590,7 @@ MODULE input_constants use_perd_yz=5,& use_perd_xyz=6,& use_perd_none=7 - + INTEGER, PARAMETER, PUBLIC :: do_clv_geo_center=0,& do_clv_fix_point =1,& do_clv_xyz =0,& @@ -598,22 +600,22 @@ MODULE input_constants do_clv_xy =4,& do_clv_xz =5,& do_clv_yz =6 - + INTEGER, PARAMETER, PUBLIC :: xc_pot_energy_none=0,& xc_pot_energy_xc_functional=1,& xc_pot_energy_sum_eigenvalues=2 - + INTEGER, PARAMETER, PUBLIC :: xc_pot_none = 1000,& xc_pot_saop = 1001 - - ! X-Ray spectra Calculation + + ! X-Ray spectra Calculation INTEGER, PARAMETER, PUBLIC :: xas_none=0,& xas_tp_hh=1,& xas_tp_fh=2,& xes_tp_val=3,& xas_tp_xhh=4,& xas_tp_xfh=5,& - xas_dscf=6 + xas_dscf=6 INTEGER, PARAMETER, PUBLIC :: xas_1s_type=1,& xas_2s_type=2,& xas_2p_type=3,& @@ -635,44 +637,44 @@ MODULE input_constants lr_chemshift=1,& lr_gtensor=2,& lr_current=3 - + INTEGER, PARAMETER, PUBLIC :: xc_pbe_orig=11,& xc_rev_pbe=12 - - ! EIP models + + ! EIP models INTEGER, PARAMETER, PUBLIC :: use_bazant_eip=1, & use_lenosky_eip=2 - - ! ddapc restraint forms + + ! ddapc restraint forms INTEGER, PARAMETER, PUBLIC :: do_ddapc_restraint=773,& do_ddapc_constraint=774 - - ! outer scf types + + ! outer scf types INTEGER, PARAMETER, PUBLIC :: outer_scf_ddapc_constraint=123, & outer_scf_s2_constraint=124,& outer_scf_becke_constraint=125, & outer_scf_scp=126, & outer_scf_none=127 - - ! outer scf optimizers + + ! outer scf optimizers INTEGER, PARAMETER, PUBLIC :: outer_scf_optimizer_sd=1001, & outer_scf_optimizer_diis=1002, & outer_scf_optimizer_none=1003, & outer_scf_optimizer_bisect=1004 - - ! s2 restraint forms + + ! s2 restraint forms INTEGER, PARAMETER, PUBLIC :: do_s2_restraint=872,& do_s2_constraint=873 - - ! ROKS schemes + + ! ROKS schemes INTEGER, PARAMETER, PUBLIC :: general_roks = 1,& high_spin_roks = 2 - - !VWN xc functional flavours + + !VWN xc functional flavours INTEGER, PARAMETER, PUBLIC :: do_vwn5=1,& do_vwn3=2 - - ! mixing of force envs + + ! mixing of force envs INTEGER, PARAMETER, PUBLIC :: mix_linear_combination=701, & mix_minimum=702, & mix_coupled=703, & @@ -686,8 +688,8 @@ MODULE input_constants ! optimization targets for distribution 2d INTEGER, PARAMETER, PUBLIC :: opt_target_minimum=456, & opt_target_median=457 - - ! HFX potential type + + ! HFX potential type INTEGER, PARAMETER, PUBLIC :: do_hfx_potential_coulomb=1,& do_hfx_potential_short=2,& do_hfx_potential_long=3,& @@ -698,11 +700,11 @@ MODULE input_constants do_hfx_potential_truncated=8 ! HFX periodic number of shells INTEGER, PARAMETER, PUBLIC :: do_hfx_auto_shells = -1 - - ! adiabatic hybrid type + + ! adiabatic hybrid type INTEGER, PARAMETER, PUBLIC :: do_adiabatic_hybrid_mcy3=1 - - ! adiabatic model type + + ! adiabatic model type INTEGER, PARAMETER, PUBLIC :: do_adiabatic_model_pade=1 ! HFX evaluation type @@ -766,7 +768,7 @@ MODULE input_constants fftw_plan_patient=3, & fftw_plan_exhaustive=4 - ! fftw will automatically be either FFTW2 or FFTW3 + ! fftw will automatically be either FFTW2 or FFTW3 #if defined ( __FFTW2 ) INTEGER, PARAMETER, PUBLIC :: do_fft_fftw=do_fft_fftw2 #else @@ -797,5 +799,5 @@ MODULE input_constants gaussian_env=2, & ramp_env=3 - + END MODULE input_constants diff --git a/src/input_cp2k_mm.F b/src/input_cp2k_mm.F index 1812ed5d87..ab36a36a9b 100644 --- a/src/input_cp2k_mm.F +++ b/src/input_cp2k_mm.F @@ -357,6 +357,10 @@ CONTAINS CALL section_add_subsection(section, subsection, error=error) CALL section_release(subsection,error=error) + CALL create_OPBEND_section(subsection, error=error) + CALL section_add_subsection(section, subsection, error=error) + CALL section_release(subsection,error=error) + CALL create_DIPOLE_section(subsection, error=error) CALL section_add_subsection(section, subsection, error=error) CALL section_release(subsection,error=error) @@ -522,7 +526,7 @@ CONTAINS END IF END SUBROUTINE create_TORSION_section - ! ***************************************************************************** +! ***************************************************************************** !> \brief This section specifies the improper torsion of the MM atoms !> \param section the section to create !> \param error variable to control error logging, stopping,... @@ -589,6 +593,75 @@ CONTAINS END IF END SUBROUTINE create_IMPROPER_section +! ***************************************************************************** +!> \brief This section specifies the out of plane bend of the MM atoms +!> \param section the section to create +!> \param error variable to control error logging, stopping,... +!> see module cp_error_handling +!> \author louis vanduyfhuys +! ***************************************************************************** + SUBROUTINE create_OPBEND_section(section,error) + TYPE(section_type), POINTER :: section + TYPE(cp_error_type), INTENT(inout) :: error + + CHARACTER(len=*), PARAMETER :: routineN = 'create_OPBEND_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="OPBEND",& + description="Specifies the out of plane bend potential of the MM system."//& + "(Only defined for atom quadruples which are also defined as an improper"//& + " pattern in the topology.)",& + n_keywords=1, n_subsections=0, repeats=.TRUE., required=.TRUE.,& + error=error) + + NULLIFY(keyword) + CALL keyword_create(keyword, name="ATOMS",& + description="Defines the atomic kinds involved in the opbend.",& + usage="ATOMS {KIND1} {KIND2} {KIND3} {KIND4}", required=.TRUE., type_of_var=char_t,& + n_var=4, error=error) + CALL section_add_keyword(section,keyword,error=error) + CALL keyword_release(keyword,error=error) + + CALL keyword_create(keyword=keyword, name="KIND",& + description="Define the kind of out of plane bend potential",& + usage="KIND HARMONIC",& + enum_c_vals=s2a("HARMONIC","MM2","MM3","MM4"),& + enum_desc=s2a("Functional Form (HARMONIC): 0.5 * K * [ PHI - PHI0 ]**2",& + "Functional Form (MM2|MM3|MM4): K * [ PHI - PHI0 ]**2",& + "Functional Form (MM2|MM3|MM4): K * [ PHI - PHI0 ]**2",& + "Functional Form (MM2|MM3|MM4): K * [ PHI - PHI0 ]**2"),& + enum_i_vals=(/ do_ff_harmonic,& + do_ff_mm2,& + do_ff_mm3,& + do_ff_mm4/),& + default_i_val=do_ff_harmonic, error=error) + CALL section_add_keyword(section,keyword,error=error) + CALL keyword_release(keyword,error=error) + + CALL keyword_create(keyword, name="K",& + description="Defines the force constant of the potential",& + usage="K {real}", required=.TRUE., type_of_var=real_t,& + n_var=1, unit_str="hartree*rad^-2",error=error) + CALL section_add_keyword(section,keyword,error=error) + CALL keyword_release(keyword,error=error) + + CALL keyword_create(keyword, name="PHI0",& + description="Defines the phase of the potential.",& + usage="PHI0 {real}", required=.TRUE., type_of_var=real_t,& + n_var=1, unit_str="rad", error=error) + CALL section_add_keyword(section,keyword,error=error) + CALL keyword_release(keyword,error=error) + + END IF + END SUBROUTINE create_OPBEND_section + ! ***************************************************************************** !> \brief This section specifies the bend of the MM atoms !> \param section the section to create diff --git a/src/mol_force.F b/src/mol_force.F index d4c4f78a62..4f41a7dd03 100644 --- a/src/mol_force.F +++ b/src/mol_force.F @@ -12,8 +12,8 @@ MODULE mol_force USE f77_blas USE input_constants, ONLY: & do_ff_amber, do_ff_charmm, do_ff_cubic, do_ff_g87, do_ff_g96, & - do_ff_harmonic, do_ff_mixed_bend_stretch, do_ff_mm3, do_ff_morse, & - do_ff_quartic + do_ff_harmonic, do_ff_mixed_bend_stretch, do_ff_mm2, do_ff_mm3, & + do_ff_mm4, do_ff_morse, do_ff_quartic USE kinds, ONLY: dp USE termination, ONLY: stop_program #include "cp_common_uses.h" @@ -22,7 +22,7 @@ MODULE mol_force PRIVATE CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mol_force' - PUBLIC :: force_bonds, force_bends, force_torsions, force_imp_torsions + PUBLIC :: force_bonds, force_bends, force_torsions, force_imp_torsions, force_opbends PUBLIC :: get_pv_bond, get_pv_bend, get_pv_torsion CONTAINS @@ -276,7 +276,7 @@ CONTAINS END SUBROUTINE force_torsions ! ***************************************************************************** -!> \brief Computes the forces from the torsions +!> \brief Computes the forces from the improper torsions !> \par History !> none !> \author DG @@ -332,6 +332,131 @@ CONTAINS gt3 = ( dist2 * is32 ** 2 - 1.0_dp ) * gt4 - dist1 * is32 ** 2 * gt1 END SUBROUTINE force_imp_torsions + ! ***************************************************************************** +!> \brief Computes the forces from the out of plane bends +!> \par History +!> none +!> \author Louis Vanduyfhuys +! ***************************************************************************** + SUBROUTINE force_opbends (id_type,s32, is32, ism, isn, tm,& + tn, t41, t42, t43, k, phi0, gt1, gt2, gt3, gt4, energy, fscalar ) + + INTEGER, INTENT(IN) :: id_type + REAL(KIND=dp), INTENT(IN) :: s32, is32, ism, isn + REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: tm, tn, t41, t42, t43 + REAL(KIND=dp), INTENT(IN) :: k, phi0 + REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: gt1, gt2, gt3, gt4 + REAL(KIND=dp), INTENT(OUT) :: energy, fscalar + + CHARACTER(len=*), PARAMETER :: routineN = 'force_opbends', & + routineP = moduleN//':'//routineN + REAL(KIND=dp), PARAMETER :: f12 = 1.0_dp/2.0_dp + + LOGICAL :: failure + REAL(KIND=dp) :: b, C, cosphi, D, E, is41, phi + REAL(KIND=dp), DIMENSION(3) :: db_dq1, db_dq2, db_dq3, & + dC_dq1, dC_dq2, dC_dq3, & + dD_dq1, dD_dq2, dD_dq3, & + dE_dq1, dE_dq2, dE_dq3 + + failure = .FALSE. + + !compute the energy and the gradients of cos(phi), see + ! "Efficient treatment of out-of-plane bend and improper torsion interactions in + ! MM2, MM3 and MM4 Molecular mechanicsd calculations.", Robert E. Tuzun, Donald W.Noid, + ! Bobby G Sumpter, Chemical and Analytical Sciences Division, Oak Ridge National + ! Laboratory, P.O. Box 2008, Oak Ridge, Tennesse 37831-6497 + + !CAUTION: in the paper r_ij = rj - ri + !in fist_intra_force.F t_ij = ri - rj + !hence a minus sign needs to be added to convert r_ij vectors in t_ij vectors + + !tm is the normal of the plane 123, tm = t32 x t12 (= w from paper) + !tn = - t41 x tm (= a from paper, for minus sign see CAUTION above) + + !Computing some necessary variables (see paper) + E = DOT_PRODUCT(-t41,tm) + C = DOT_PRODUCT(tm,tm) + D = E**2/C + b = DOT_PRODUCT(t41,t41)-D + + !inverse norm of t41 + is41 = 1.0_dp/SQRT(DOT_PRODUCT(t41,t41)) + + cosphi = SQRT(b)*is41 + IF(cosphi > 1.0_dp) cosphi = 1.0_dp + IF(cosphi < -1.0_dp) cosphi = -1.0_dp + phi = SIGN( ACOS(cosphi) , DOT_PRODUCT(tm,-t41) ) + + SELECT CASE (id_type) + CASE (do_ff_mm2,do_ff_mm3,do_ff_mm4) + ! compute energy + energy = k * ( phi - phi0 )**2 + + ! compute fscalar + fscalar = 2.0_dp * k * ( phi - phi0 ) * is41 + + CASE (do_ff_harmonic) + ! compute energy + energy = f12 * k * ( phi - phi0 )**2 + + ! compute fscalar + fscalar = k * ( phi - phi0 ) * is41 + + CASE DEFAULT + CALL stop_program (routineP,"unmatched opbend kind") + END SELECT + + + !Computing the necessary intermediate variables. dX_dqi is the gradient + !of X with respect to the coördinates of particle i. + + dE_dq1(1) = ( t42(2)*t43(3)-t43(2)*t42(3)) + dE_dq1(2) = (-t42(1)*t43(3)+t43(1)*t42(3)) + dE_dq1(3) = ( t42(1)*t43(2)-t43(1)*t42(2)) + + dE_dq2(1) = ( t43(2)*t41(3)-t41(2)*t43(3)) + dE_dq2(2) = (-t43(1)*t41(3)+t41(1)*t43(3)) + dE_dq2(3) = ( t43(1)*t41(2)-t41(1)*t43(2)) + + dE_dq3(1) = ( t41(2)*t42(3)-t42(2)*t41(3)) + dE_dq3(2) = (-t41(1)*t42(3)+t42(1)*t41(3)) + dE_dq3(3) = ( t41(1)*t42(2)-t42(1)*t41(2)) + + + dC_dq1 = 2.0_dp * ( (t42-t41)*s32**2 - (t42-t43)*DOT_PRODUCT(t42-t41,t42-t43) ) + dC_dq3 = 2.0_dp * ( (t42-t43)*DOT_PRODUCT(t41-t42,t41-t42) & + -(t42-t41)*DOT_PRODUCT(t42-t41,t42-t43) ) + !C only dependent of atom 1 2 and 3, using translational invariance we find + dC_dq2 = - ( dC_dq1 + dC_dq3 ) + + dD_dq1 = (2.0_dp*E*dE_dq1 - D*dC_dq1)/C + dD_dq2 = (2.0_dp*E*dE_dq2 - D*dC_dq2)/C + dD_dq3 = (2.0_dp*E*dE_dq3 - D*dC_dq3)/C + + db_dq1 = - 2.0_dp*t41 - dD_dq1 + db_dq2 = - dD_dq2 + db_dq3 = - dD_dq3 + + !gradients of cos(phi), gt4 is calculated using translational invariance. + !The 1/r41 factor from the paper is absorbed in fscalar. + !If phi = 0 then sin(phi) = 0 and the regular formula for calculating gt + !won't work because of the sine function in the denominator. A further + !analytic simplification is needed. + IF (phi == 0) THEN + gt1 = - SIGN(1.0_dp,phi)/SQRT(C)*dE_dq1 + gt2 = - SIGN(1.0_dp,phi)/SQRT(C)*dE_dq2 + gt3 = - SIGN(1.0_dp,phi)/SQRT(C)*dE_dq3 + gt4 = - (gt1 + gt2 + gt3) + + ELSE + gt1 = (1.0_dp/(2.0_dp*SQRT(b))*db_dq1 + cosphi*t41*is41)/SIN(phi) + gt2 = (1.0_dp/(2.0_dp*SQRT(b))*db_dq2)/SIN(phi) + gt3 = (1.0_dp/(2.0_dp*SQRT(b))*db_dq3)/SIN(phi) + gt4 = - (gt1 + gt2 + gt3) + END IF + END SUBROUTINE force_opbends + ! ***************************************************************************** !> \brief Computes the pressure tensor from the bonds !> \par History @@ -394,7 +519,8 @@ CONTAINS END SUBROUTINE get_pv_bend ! ***************************************************************************** -!> \brief Computes the pressure tensor from the torsions +!> \brief Computes the pressure tensor from the torsions (also used for impr +!> and opbend) !> \par History !> none !> \author DG @@ -437,5 +563,6 @@ CONTAINS END SUBROUTINE get_pv_torsion + END MODULE mol_force diff --git a/src/molecule_kind_types.F b/src/molecule_kind_types.F index 80d166c6a8..d3f1801134 100644 --- a/src/molecule_kind_types.F +++ b/src/molecule_kind_types.F @@ -30,8 +30,9 @@ MODULE molecule_kind_types USE f77_blas USE force_field_types, ONLY: & bend_kind_type, bond_kind_type, impr_kind_dealloc_ref, impr_kind_type, & - onfo_kind_dealloc_ref, onfo_kind_type, torsion_kind_dealloc_ref, & - torsion_kind_type, ub_kind_dealloc_ref, ub_kind_type + onfo_kind_dealloc_ref, onfo_kind_type, opbend_kind_type, & + torsion_kind_dealloc_ref, torsion_kind_type, ub_kind_dealloc_ref, & + ub_kind_type USE input_constants, ONLY: use_perd_x,& use_perd_xy,& use_perd_xyz,& @@ -102,6 +103,13 @@ MODULE molecule_kind_types TYPE(impr_kind_type), POINTER :: impr_kind END TYPE impr_type +! ***************************************************************************** + TYPE opbend_type + INTEGER :: a,b,c,d + INTEGER :: id_type, itype + TYPE(opbend_kind_type), POINTER :: opbend_kind + END TYPE opbend_type + ! ***************************************************************************** TYPE onfo_type INTEGER :: a,b @@ -172,6 +180,8 @@ MODULE molecule_kind_types TYPE(torsion_type), DIMENSION(:), POINTER :: torsion_list TYPE(impr_kind_type), DIMENSION(:), POINTER :: impr_kind_set TYPE(impr_type), DIMENSION(:), POINTER :: impr_list + TYPE(opbend_kind_type), DIMENSION(:), POINTER :: opbend_kind_set + TYPE(opbend_type), DIMENSION(:), POINTER :: opbend_list TYPE(onfo_kind_type), DIMENSION(:), POINTER :: onfo_kind_set TYPE(onfo_type), DIMENSION(:), POINTER :: onfo_list TYPE(colvar_constraint_type), DIMENSION(:), POINTER :: colv_list @@ -188,6 +198,7 @@ MODULE molecule_kind_types nbond,& nbend,& nimpr,& + nopbend,& nonfo,& ntorsion,& nub,& @@ -220,6 +231,7 @@ MODULE molecule_kind_types ub_type,& torsion_type,& impr_type,& + opbend_type,& onfo_type,& restraint_type,& colvar_constraint_type,& @@ -353,6 +365,8 @@ CONTAINS NULLIFY (molecule_kind_set(imolecule_kind)%ub_kind_set) NULLIFY (molecule_kind_set(imolecule_kind)%impr_kind_set) NULLIFY (molecule_kind_set(imolecule_kind)%impr_list) + NULLIFY (molecule_kind_set(imolecule_kind)%opbend_kind_set) + NULLIFY (molecule_kind_set(imolecule_kind)%opbend_list) NULLIFY (molecule_kind_set(imolecule_kind)%onfo_list) NULLIFY (molecule_kind_set(imolecule_kind)%onfo_kind_set) NULLIFY (molecule_kind_set(imolecule_kind)%g3x3_list) @@ -374,6 +388,7 @@ CONTAINS molecule_kind_set(imolecule_kind)%nbond = 0 molecule_kind_set(imolecule_kind)%nimpr = 0 molecule_kind_set(imolecule_kind)%nonfo = 0 + molecule_kind_set(imolecule_kind)%nopbend = 0 molecule_kind_set(imolecule_kind)%nub = 0 CALL setup_colvar_counters(molecule_kind_set(imolecule_kind)%colv_list,& molecule_kind_set(imolecule_kind)%ncolv,error) @@ -449,6 +464,14 @@ CONTAINS DEALLOCATE (molecule_kind_set(imolecule_kind)%impr_kind_set,STAT=stat) CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) END IF + IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%opbend_list)) THEN + DEALLOCATE (molecule_kind_set(imolecule_kind)%opbend_list,STAT=stat) + CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) + END IF + IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%opbend_kind_set)) THEN + DEALLOCATE (molecule_kind_set(imolecule_kind)%opbend_kind_set,STAT=stat) + CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) + END IF IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%onfo_list)) THEN DEALLOCATE(molecule_kind_set(imolecule_kind)%onfo_list, STAT=stat) CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) @@ -525,14 +548,14 @@ CONTAINS !> \version 1.0 ! ***************************************************************************** SUBROUTINE get_molecule_kind(molecule_kind,atom_list,bond_list,bend_list,& - ub_list,impr_list,onfo_list,colv_list,fixd_list,& + ub_list,impr_list,opbend_list,onfo_list,colv_list,fixd_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,& + nimpr,nopbend,nonfo,nconstraint,nconstraint_fixd,nfixd,ncolv,ng3x3,ng4x6,& 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,& + ub_kind_set,impr_kind_set,opbend_kind_set,onfo_kind_set,torsion_kind_set,& molname_generated) TYPE(molecule_kind_type), POINTER :: molecule_kind @@ -546,6 +569,8 @@ CONTAINS POINTER :: ub_list TYPE(impr_type), DIMENSION(:), & OPTIONAL, POINTER :: impr_list + TYPE(opbend_type), DIMENSION(:), & + OPTIONAL, POINTER :: opbend_list TYPE(onfo_type), DIMENSION(:), & OPTIONAL, POINTER :: onfo_list TYPE(colvar_constraint_type), & @@ -565,10 +590,8 @@ CONTAINS CHARACTER(LEN=default_string_length), & INTENT(OUT), OPTIONAL :: name REAL(KIND=dp), OPTIONAL :: mass, charge - INTEGER, INTENT(OUT), OPTIONAL :: kind_number, natom, nbend, & - nbond, nub, nimpr, nonfo, & - nconstraint, & - nconstraint_fixd, nfixd + INTEGER, INTENT(OUT), OPTIONAL :: kind_number, natom, nbend, nbond, nub, & + nimpr, nopbend, nonfo, nconstraint, nconstraint_fixd, nfixd TYPE(colvar_counters), INTENT(out), & OPTIONAL :: ncolv INTEGER, INTENT(OUT), OPTIONAL :: ng3x3, ng4x6, nvsite, nfixd_restraint, & @@ -584,6 +607,8 @@ CONTAINS OPTIONAL, POINTER :: ub_kind_set TYPE(impr_kind_type), DIMENSION(:), & OPTIONAL, POINTER :: impr_kind_set + TYPE(opbend_kind_type), DIMENSION(:), & + OPTIONAL, POINTER :: opbend_kind_set TYPE(onfo_kind_type), DIMENSION(:), & OPTIONAL, POINTER :: onfo_kind_set TYPE(torsion_kind_type), DIMENSION(:), & @@ -601,12 +626,14 @@ CONTAINS IF (PRESENT(bend_list)) bend_list => molecule_kind%bend_list IF (PRESENT(bond_list)) bond_list => molecule_kind%bond_list IF (PRESENT(impr_list)) impr_list => molecule_kind%impr_list + IF (PRESENT(opbend_list)) opbend_list => molecule_kind%opbend_list IF (PRESENT(onfo_list)) onfo_list => molecule_kind%onfo_list IF (PRESENT(ub_list)) ub_list => molecule_kind%ub_list IF (PRESENT(bond_kind_set)) bond_kind_set => molecule_kind%bond_kind_set IF (PRESENT(bend_kind_set)) bend_kind_set => molecule_kind%bend_kind_set IF (PRESENT(ub_kind_set)) ub_kind_set => molecule_kind%ub_kind_set IF (PRESENT(impr_kind_set)) impr_kind_set => molecule_kind%impr_kind_set + IF (PRESENT(opbend_kind_set)) opbend_kind_set => molecule_kind%opbend_kind_set IF (PRESENT(onfo_kind_set)) onfo_kind_set => molecule_kind%onfo_kind_set IF (PRESENT(torsion_kind_set)) torsion_kind_set => molecule_kind%torsion_kind_set IF (PRESENT(colv_list)) colv_list => molecule_kind%colv_list @@ -626,6 +653,7 @@ CONTAINS IF (PRESENT(nbond)) nbond = molecule_kind%nbond IF (PRESENT(nub)) nub = molecule_kind%nub IF (PRESENT(nimpr)) nimpr = molecule_kind%nimpr + IF (PRESENT(nopbend)) nopbend = molecule_kind%nopbend 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) +& @@ -638,7 +666,7 @@ CONTAINS ! Number of atoms that have one or more components fixed IF (PRESENT(nfixd)) nfixd= molecule_kind%nfixd ! Number of degrees of freedom fixed - IF (PRESENT(nconstraint_fixd)) THEN + IF (PRESENT(nconstraint_fixd)) THEN nconstraint_fixd = 0 IF (molecule_kind%nfixd/=0) THEN DO i = 1, SIZE(molecule_kind%fixd_list) @@ -685,22 +713,20 @@ CONTAINS !> \version 1.0 ! ***************************************************************************** SUBROUTINE get_molecule_kind_set(molecule_kind_set,maxatom,natom,& - nbond,nbend,nub,ntorsion,nimpr,nonfo,& + nbond,nbend,nub,ntorsion,nimpr,nopbend,nonfo,& nconstraint,nmolecule,nrestraints) TYPE(molecule_kind_type), DIMENSION(:), & POINTER :: molecule_kind_set - INTEGER, INTENT(OUT), OPTIONAL :: maxatom, natom, nbond, nbend, & - nub, ntorsion, nimpr, nonfo, & - nconstraint, nmolecule, & - nrestraints + INTEGER, INTENT(OUT), OPTIONAL :: maxatom, natom, nbond, nbend, nub, & + ntorsion, nimpr, nopbend, nonfo, nconstraint, nmolecule, nrestraints CHARACTER(len=*), PARAMETER :: routineN = 'get_molecule_kind_set', & routineP = moduleN//':'//routineN - INTEGER :: ibend, ibond, iimpr, imolecule_kind, ionfo, itorsion, iub, na, & - nc, nconstraint_fixd, nfixd, nfixd_restraint, nm, nmolecule_kind, & - nrestraints_tot + INTEGER :: ibend, ibond, iimpr, imolecule_kind, ionfo, iopbend, itorsion, & + iub, na, nc, nconstraint_fixd, nfixd, nfixd_restraint, nm, & + nmolecule_kind, nrestraints_tot TYPE(molecule_kind_type), POINTER :: molecule_kind IF (ASSOCIATED(molecule_kind_set)) THEN @@ -712,7 +738,7 @@ CONTAINS IF (PRESENT(nub)) nub = 0 IF (PRESENT(ntorsion)) ntorsion = 0 IF (PRESENT(nimpr)) nimpr = 0 - IF (PRESENT(nimpr)) nimpr = 0 + IF (PRESENT(nopbend)) nopbend = 0 IF (PRESENT(nonfo)) nonfo = 0 IF (PRESENT(nconstraint)) nconstraint = 0 IF (PRESENT(nrestraints)) nrestraints = 0 @@ -731,6 +757,7 @@ CONTAINS nub=iub,& ntorsion=itorsion,& nimpr=iimpr,& + nopbend=iopbend,& nonfo=ionfo,& nconstraint=nc,& nconstraint_fixd=nconstraint_fixd,& @@ -745,6 +772,7 @@ CONTAINS IF (PRESENT(nub)) nub = nub + iub*nm IF (PRESENT(ntorsion)) ntorsion = ntorsion + itorsion*nm IF (PRESENT(nimpr)) nimpr = nimpr + iimpr*nm + IF (PRESENT(nopbend)) nopbend = nopbend + iopbend*nm IF (PRESENT(nonfo)) nonfo = nonfo + ionfo*nm IF (PRESENT(nconstraint)) nconstraint = nconstraint + nc*nm + nconstraint_fixd IF (PRESENT(nmolecule)) nmolecule = nmolecule + nm @@ -844,13 +872,15 @@ CONTAINS SUBROUTINE set_molecule_kind(molecule_kind,name,mass,charge,kind_number,& molecule_list,atom_list,nbond,bond_list,& 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,& - 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) + nopbend, opbend_list, nonfo,onfo_list,ntorsion,& + torsion_list,fixd_list,ncolv,colv_list,ng3x3,& + g3x3_list,ng4x6,nfixd,g4x6_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,& + opbend_kind_set,onfo_kind_set,nelectron,nsgf,& + molname_generated) TYPE(molecule_kind_type), POINTER :: molecule_kind CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: name @@ -871,6 +901,9 @@ CONTAINS INTEGER, INTENT(IN), OPTIONAL :: nimpr TYPE(impr_type), DIMENSION(:), & OPTIONAL, POINTER :: impr_list + INTEGER, INTENT(IN), OPTIONAL :: nopbend + TYPE(opbend_type), DIMENSION(:), & + OPTIONAL, POINTER :: opbend_list INTEGER, INTENT(IN), OPTIONAL :: nonfo TYPE(onfo_type), DIMENSION(:), & OPTIONAL, POINTER :: onfo_list @@ -908,6 +941,8 @@ CONTAINS OPTIONAL, POINTER :: torsion_kind_set TYPE(impr_kind_type), DIMENSION(:), & OPTIONAL, POINTER :: impr_kind_set + TYPE(opbend_kind_type), DIMENSION(:), & + OPTIONAL, POINTER :: opbend_kind_set TYPE(onfo_kind_type), DIMENSION(:), & OPTIONAL, POINTER :: onfo_kind_set INTEGER, INTENT(IN), OPTIONAL :: nelectron, nsgf @@ -942,6 +977,8 @@ CONTAINS IF (PRESENT(torsion_list)) molecule_kind%torsion_list => torsion_list IF (PRESENT(nimpr)) molecule_kind%nimpr = nimpr IF (PRESENT(impr_list)) molecule_kind%impr_list => impr_list + IF (PRESENT(nopbend)) molecule_kind%nopbend = nopbend + IF (PRESENT(opbend_list)) molecule_kind%opbend_list => opbend_list IF (PRESENT(nonfo)) molecule_kind%nonfo = nonfo IF (PRESENT(onfo_list)) molecule_kind%onfo_list => onfo_list IF (PRESENT(ncolv)) molecule_kind%ncolv = ncolv @@ -963,6 +1000,7 @@ CONTAINS IF (PRESENT(ub_kind_set)) molecule_kind%ub_kind_set => ub_kind_set IF (PRESENT(torsion_kind_set)) molecule_kind%torsion_kind_set => torsion_kind_set IF (PRESENT(impr_kind_set)) molecule_kind%impr_kind_set => impr_kind_set + IF (PRESENT(opbend_kind_set)) molecule_kind%opbend_kind_set => opbend_kind_set IF (PRESENT(onfo_kind_set)) molecule_kind%onfo_kind_set => onfo_kind_set IF (PRESENT(nshell)) molecule_kind%nshell = nshell IF (PRESENT(shell_list)) molecule_kind%shell_list => shell_list @@ -1051,6 +1089,9 @@ CONTAINS IF ( molecule_kind%nimpr > 0 ) & WRITE (UNIT=output_unit,FMT="(1X,A30,I6)")& "Number of improper: ",molecule_kind%nimpr + IF ( molecule_kind%nopbend > 0 ) & + WRITE (UNIT=output_unit,FMT="(1X,A30,I6)")& + "Number of out opbends: ",molecule_kind%nopbend IF ( molecule_kind%nonfo > 0 ) & WRITE (UNIT=output_unit,FMT="(1X,A30,I6)")& "Number of one-four: ",molecule_kind%nonfo @@ -1078,10 +1119,8 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'write_molecule_kind_set', & routineP = moduleN//':'//routineN - INTEGER :: imolecule_kind, nbend, nbond, & - nimpr, nmolecule_kind, nonfo, & - ntors, ntotal, nub, & - output_unit + INTEGER :: imolecule_kind, nbend, nbond, nimpr, nmolecule_kind, nonfo, & + nopbend, ntors, ntotal, nub, output_unit LOGICAL :: failure TYPE(cp_logger_type), POINTER :: logger TYPE(molecule_kind_type), POINTER :: molecule_kind @@ -1108,8 +1147,9 @@ CONTAINS nub=nub,& ntorsion=ntors,& nimpr=nimpr,& + nopbend=nopbend,& nonfo=nonfo) - ntotal = nbond+nbend+nub+ntors+nimpr+nonfo + ntotal = nbond+nbend+nub+ntors+nimpr+nopbend+nonfo IF ( ntotal > 0 ) THEN WRITE(UNIT=output_unit,FMT="(/,/,T2,A,T45,A30,I6)") & "MOLECULE KIND SET INFORMATION", & @@ -1122,6 +1162,8 @@ CONTAINS "Total Number of torsions: ",ntors WRITE (UNIT=output_unit,FMT="(T45,A30,I6)")& "Total Number of improper: ",nimpr + WRITE (UNIT=output_unit,FMT="(T45,A30,I6)")& + "Total Number of opbends: ",nopbend WRITE (UNIT=output_unit,FMT="(T45,A30,I6)")& "Total Number of one-four: ",nonfo END IF diff --git a/src/topology_connectivity_util.F b/src/topology_connectivity_util.F index 7792b9749d..da6b5aa8a5 100644 --- a/src/topology_connectivity_util.F +++ b/src/topology_connectivity_util.F @@ -8,13 +8,15 @@ !> \par History !> jgh (23-05-2004) Last atom of molecule information added ! ***************************************************************************** + MODULE topology_connectivity_util USE cp_output_handling, ONLY: cp_print_key_finished_output,& cp_print_key_unit_nr USE f77_blas USE input_constants, ONLY: do_conn_g87,& do_conn_g96,& - do_ff_charmm + do_ff_charmm,& + do_ff_harmonic USE input_section_types, ONLY: section_vals_type,& section_vals_val_get USE kinds, ONLY: default_string_length @@ -22,7 +24,7 @@ MODULE topology_connectivity_util USE molecule_kind_types, ONLY: & allocate_molecule_kind_set, atom_type, bend_type, bond_type, & get_molecule_kind, impr_type, molecule_kind_type, onfo_type, & - set_molecule_kind, torsion_type, ub_type + opbend_type, set_molecule_kind, torsion_type, ub_type USE molecule_types_new, ONLY: allocate_molecule_set,& get_molecule,& local_molecule_type,& @@ -49,6 +51,8 @@ CONTAINS ! ***************************************************************************** !> \brief topology connectivity pack +!> \History 11/2009 (Louis Vanduyhuys): added Out of Plane bends based on +!> impropers in topology ! ***************************************************************************** SUBROUTINE topology_connectivity_pack(molecule_kind_set,molecule_set,& topology,subsys_section,error) @@ -88,6 +92,7 @@ CONTAINS TYPE(molecule_kind_type), POINTER :: molecule_kind TYPE(molecule_type), POINTER :: molecule TYPE(onfo_type), DIMENSION(:), POINTER :: onfo_list + TYPE(opbend_type), DIMENSION(:), POINTER :: opbend_list TYPE(torsion_type), DIMENSION(:), & POINTER :: torsion_list TYPE(ub_type), DIMENSION(:), POINTER :: ub_list @@ -100,7 +105,7 @@ CONTAINS output_unit= cp_logger_get_default_unit_nr(logger) ELSE output_unit = -1 - END IF + END IF iw = cp_print_key_unit_nr(logger,subsys_section,"PRINT%TOPOLOGY_INFO/UTIL_INFO",& extension=".subsysLog",error=error) CALL timeset(routineN,handle) @@ -410,7 +415,7 @@ CONTAINS END IF molecule_kind => molecule_kind_set(i) nval_tot2 = nval_tot2 + ibond*SIZE(molecule_kind%molecule_list) - + ALLOCATE(bond_list(ibond),STAT=stat) CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) ibond = 0 @@ -1030,6 +1035,7 @@ CONTAINS !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- ! 14. Set the molecule_kind%[nimpr,impr_list] via set_molecule_kind + ! Also set the molecule_kind%[nopbend,opbend_list] !----------------------------------------------------------------------------- ! Allocate c_var_a, c_var_b, c_var_c, c_var_d, c_var_type CALL timeset(routineN//"_14_pre",handle2) @@ -1115,41 +1121,56 @@ CONTAINS CALL find_bnd_typ(topology%nmol_type, map_vars, map_var_mol, bnd_type, nvar1, error) CALL find_bnd_typ(topology%nmol_type, map_cvars, map_cvar_mol, bnd_ctype, nvar2, error) DO i=1,topology%nmol_type - intra_imprs = 0 - inter_imprs = 0 + intra_imprs = 0 + inter_imprs = 0 IF (ALL(bnd_type(:,i)>0)) THEN - intra_imprs = bnd_type(2,i)-bnd_type(1,i)+1 + intra_imprs = bnd_type(2,i)-bnd_type(1,i)+1 END IF IF (ALL(bnd_ctype(:,i)>0)) THEN - inter_imprs = bnd_ctype(2,i)-bnd_ctype(1,i)+1 + inter_imprs = bnd_ctype(2,i)-bnd_ctype(1,i)+1 END IF - iimpr = intra_imprs+inter_imprs + iimpr = intra_imprs+inter_imprs IF(iw>0) THEN WRITE(iw,*) " Total number of imprs for molecule type ",i," :",iimpr WRITE(iw,*) " intra (imprs inside molecules) :: ",intra_imprs WRITE(iw,*) " inter (imprs between molecules) :: ",inter_imprs + WRITE(iw,*) " Total number of opbends for molecule type ",i," :",iimpr + WRITE(iw,*) " intra (opbends inside molecules) :: ",intra_imprs + WRITE(iw,*) " inter (opbends between molecules) :: ",inter_imprs END IF molecule_kind => molecule_kind_set(i) nval_tot2 = nval_tot2 + iimpr*SIZE(molecule_kind%molecule_list) ALLOCATE(impr_list(iimpr),STAT=stat) + ALLOCATE(opbend_list(iimpr),STAT=stat) CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) - iimpr = 0 + iimpr = 0 DO j=bnd_type(1,i),bnd_type(2,i) IF (j==0) CYCLE - iimpr = iimpr + 1 + iimpr = iimpr + 1 jind = map_vars(j) first = first_list(map_atom_mol(conn_info%impr_a(jind))) impr_list(iimpr)%a=conn_info%impr_a(jind)-first+1 impr_list(iimpr)%b=conn_info%impr_b(jind)-first+1 impr_list(iimpr)%c=conn_info%impr_c(jind)-first+1 impr_list(iimpr)%d=conn_info%impr_d(jind)-first+1 - ! Set by default id_type to charmm and modify when handling the forcefield + ! Atom sequence for improper is A B C D in which A is central atom, + ! B is deviating atom and C & D are secondairy atoms. Atom sequence for + ! opbend is B D C A in which A is central atom, B is deviating. Hence + ! to create an opbend out of an improper, B and D need to be interchanged. + opbend_list(iimpr)%a=conn_info%impr_b(jind)-first+1 + opbend_list(iimpr)%b=conn_info%impr_d(jind)-first+1 + opbend_list(iimpr)%c=conn_info%impr_c(jind)-first+1 + opbend_list(iimpr)%d=conn_info%impr_a(jind)-first+1 + ! Set by default id_type of improper to charmm and modify when handling the forcefield impr_list(iimpr)%id_type = do_ff_charmm + ! Set by default id_type of opbend to harmonic and modify when handling the forcefield + opbend_list(iimpr)%id_type = do_ff_harmonic IF((topology%conn_type==do_conn_g96).OR.(topology%conn_type==do_conn_g87)) THEN impr_list(iimpr)%itype = conn_info%impr_type(jind) END IF !point this to the right impr_kind_type if using force field NULLIFY(impr_list(iimpr)%impr_kind) + NULLIFY(opbend_list(iimpr)%opbend_kind) IF(iw>0) THEN WRITE(iw,'(7X,A,I3,1X,A,I4,I4,I4,I4,1X,A,I4,I4,I4,I4)') & "molecule_kind",i,"intra IMPR",& @@ -1162,6 +1183,17 @@ CONTAINS conn_info%impr_b(jind)-first+1,& conn_info%impr_c(jind)-first+1,& conn_info%impr_d(jind)-first+1 + WRITE(iw,'(7X,A,I3,1X,A,I4,I4,I4,I4,1X,A,I4,I4,I4,I4)') & + "molecule_kind",i,"intra OPBEND",& + conn_info%impr_b(jind),& + conn_info%impr_d(jind),& + conn_info%impr_c(jind),& + conn_info%impr_a(jind),& + "offset number at",& + conn_info%impr_b(jind)-first+1,& + conn_info%impr_d(jind)-first+1,& + conn_info%impr_c(jind)-first+1,& + conn_info%impr_a(jind)-first+1 END IF END DO DO j=bnd_ctype(1,i),bnd_ctype(2,i) @@ -1174,13 +1206,20 @@ CONTAINS impr_list(iimpr)%b=c_var_b(jind)-first+1 impr_list(iimpr)%c=c_var_c(jind)-first+1 impr_list(iimpr)%d=c_var_d(jind)-first+1 - ! Set by default id_type to charmm and modify when handling the forcefield + opbend_list(iimpr)%a=c_var_b(jind)-first+1 + opbend_list(iimpr)%b=c_var_d(jind)-first+1 + opbend_list(iimpr)%c=c_var_c(jind)-first+1 + opbend_list(iimpr)%d=c_var_a(jind)-first+1 + ! Set by default id_type of improper to charmm and modify when handling the forcefield impr_list(iimpr)%id_type = do_ff_charmm + ! Set by default id_type of opbend to harmonic and modify when handling the forcefield + opbend_list(iimpr)%id_type = do_ff_harmonic IF((topology%conn_type==do_conn_g96).OR.(topology%conn_type==do_conn_g87)) THEN impr_list(iimpr)%itype = c_var_type(jind) END IF - !point this to the right impr_kind_type if using force field + !point this to the right impr_kind_type and opbend_kind_type if using force field NULLIFY(impr_list(iimpr)%impr_kind) + NULLIFY(opbend_list(iimpr)%opbend_kind) IF(iw>0) THEN WRITE(iw,'(7X,A,I3,1X,A,I4,I4,I4,I4,1X,A,I4,I4,I4,I4)') & "molecule_kind",i,"inter IMPR",& @@ -1193,10 +1232,23 @@ CONTAINS c_var_b(jind)-first+1,& c_var_c(jind)-first+1,& c_var_d(jind)-first+1 + WRITE(iw,'(7X,A,I3,1X,A,I4,I4,I4,I4,1X,A,I4,I4,I4,I4)') & + "molecule_kind",i,"inter OPBEND",& + c_var_b(jind),& + c_var_d(jind),& + c_var_c(jind),& + c_var_a(jind),& + "offset number at",& + c_var_b(jind)-first+1,& + c_var_d(jind)-first+1,& + c_var_c(jind)-first+1,& + c_var_a(jind)-first+1 END IF END DO CALL set_molecule_kind(molecule_kind=molecule_kind,& nimpr=SIZE(impr_list),impr_list=impr_list) + CALL set_molecule_kind(molecule_kind=molecule_kind,& + nopbend=SIZE(opbend_list),opbend_list=opbend_list) END DO CPPostcondition(nval_tot1==nval_tot2,cp_failure_level,routineP,error,failure) DEALLOCATE(map_var_mol,stat=stat) diff --git a/tests/Fist/regtest-14/TEST_FILES b/tests/Fist/regtest-14/TEST_FILES index 145eb375ec..a2871c246a 100644 --- a/tests/Fist/regtest-14/TEST_FILES +++ b/tests/Fist/regtest-14/TEST_FILES @@ -13,3 +13,4 @@ imp_test_11.inp 11 imp_test_12.inp 11 imp_test_21.inp 11 imp_test_22.inp 11 +opbend_test.inp 11 diff --git a/tests/Fist/regtest-14/imp_test_12.inp b/tests/Fist/regtest-14/imp_test_12.inp index e94e248a1a..6ec0d28c35 100644 --- a/tests/Fist/regtest-14/imp_test_12.inp +++ b/tests/Fist/regtest-14/imp_test_12.inp @@ -178,9 +178,8 @@ H 0.10430 0.07370 1.01100 &END COORD &TOPOLOGY - &GENERATE - CREATE_MOLECULES - &END GENERATE + CONNECTIVITY UPSF + CONN_FILE_NAME ../sample_psf/imp_test.psf &END TOPOLOGY &END SUBSYS &END FORCE_EVAL diff --git a/tests/Fist/regtest-14/opbend_test.inp b/tests/Fist/regtest-14/opbend_test.inp new file mode 100644 index 0000000000..5b2264041a --- /dev/null +++ b/tests/Fist/regtest-14/opbend_test.inp @@ -0,0 +1,178 @@ +&FORCE_EVAL + METHOD FIST + &MM + &FORCEFIELD + &BOND + ATOMS C H + K [au_e*au_l^-2] 0.0 + R0 [au_l] 1.91807202593 + &END BOND + &BOND + ATOMS C Cl + K [au_e*au_l^-2] 0.0 + R0 [au_l] 1.91807202593 + &END BOND + &BOND + ATOMS C Br + K [au_e*au_l^-2] 0.0 + R0 [au_l] 1.91807202593 + &END BOND + &BEND + ATOMS Br C H + K [au_e*rad^-2] 0.0 + THETA0 [rad] 0.0 + &END BEND + &BEND + ATOMS Br C Cl + K [au_e*rad^-2] 0.0 + THETA0 [rad] 0.0 + &END BEND + &BEND + ATOMS Cl C H + K [au_e*rad^-2] 0.0 + THETA0 [rad] 0.0 + &END BEND + &OPBEND + KIND MM3 + ATOMS Br H Cl C + K [au_e*rad^-2] 2.0 + PHI0 [rad] 1.1123 + &END OPBEND + &NONBONDED + &LENNARD-JONES + ATOMS C C + SIGMA [au_l] 1.0 + EPSILON [au_e] 0.0 + &END LENNARD-JONES + &LENNARD-JONES + ATOMS C H + SIGMA [au_l] 1.0 + EPSILON [au_e] 0.0 + &END LENNARD-JONES + &LENNARD-JONES + ATOMS C Cl + SIGMA [au_l] 1.0 + EPSILON [au_e] 0.0 + &END LENNARD-JONES + &LENNARD-JONES + ATOMS C Br + SIGMA [au_l] 1.0 + EPSILON [au_e] 0.0 + &END LENNARD-JONES + &LENNARD-JONES + ATOMS H H + SIGMA [au_l] 1.0 + EPSILON [au_e] 0.0 + &END LENNARD-JONES + &LENNARD-JONES + ATOMS H Br + SIGMA [au_l] 1.0 + EPSILON [au_e] 0.0 + &END LENNARD-JONES + &LENNARD-JONES + ATOMS H Cl + SIGMA [au_l] 1.0 + EPSILON [au_e] 0.0 + &END LENNARD-JONES + &LENNARD-JONES + ATOMS Br Br + SIGMA [au_l] 1.0 + EPSILON [au_e] 0.0 + &END LENNARD-JONES + &LENNARD-JONES + ATOMS Br Cl + SIGMA [au_l] 1.0 + EPSILON [au_e] 0.0 + &END LENNARD-JONES + &LENNARD-JONES + ATOMS Cl Cl + SIGMA [au_l] 1.0 + EPSILON [au_e] 0.0 + &END LENNARD-JONES + &END NONBONDED + &NONBONDED14 + &LENNARD-JONES + ATOMS C C + SIGMA [au_l] 1.0 + EPSILON [au_e] 0.0 + &END LENNARD-JONES + &LENNARD-JONES + ATOMS C H + SIGMA [au_l] 1.0 + EPSILON [au_e] 0.0 + &END LENNARD-JONES + &LENNARD-JONES + ATOMS C Cl + SIGMA [au_l] 1.0 + EPSILON [au_e] 0.0 + &END LENNARD-JONES + &LENNARD-JONES + ATOMS C Br + SIGMA [au_l] 1.0 + EPSILON [au_e] 0.0 + &END LENNARD-JONES + &LENNARD-JONES + ATOMS H H + SIGMA [au_l] 1.0 + EPSILON [au_e] 0.0 + &END LENNARD-JONES + &LENNARD-JONES + ATOMS H Br + SIGMA [au_l] 1.0 + EPSILON [au_e] 0.0 + &END LENNARD-JONES + &LENNARD-JONES + ATOMS H Cl + SIGMA [au_l] 1.0 + EPSILON [au_e] 0.0 + &END LENNARD-JONES + &LENNARD-JONES + ATOMS Br Br + SIGMA [au_l] 1.0 + EPSILON [au_e] 0.0 + &END LENNARD-JONES + &LENNARD-JONES + ATOMS Br Cl + SIGMA [au_l] 1.0 + EPSILON [au_e] 0.0 + &END LENNARD-JONES + &LENNARD-JONES + ATOMS Cl Cl + SIGMA [au_l] 1.0 + EPSILON [au_e] 0.0 + &END LENNARD-JONES + &END NONBONDED14 + &END FORCEFIELD + &POISSON + PERIODIC none + &EWALD + EWALD_TYPE none + &END EWALD + &END POISSON + &END MM + &SUBSYS + &CELL + ABC [angstrom] 25.9457 25.9457 25.9457 + PERIODIC NONE + &END CELL + &COORD + Br -1.45470 0.99560 -0.57390 + Cl -0.18290 -1.68550 -0.47690 + C 0.00000 0.00000 0.00000 + H 0.10430 0.07370 1.01100 + &END COORD + &TOPOLOGY + CONNECTIVITY UPSF + CONN_FILE_NAME ../sample_psf/bromochloroamine.psf + &END TOPOLOGY + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT opbend_test + RUN_TYPE DEBUG + PRINT_LEVEL DEBUG +&END GLOBAL +&DEBUG + DX 0.001 + DEBUG_STRESS_TENSOR F +&END diff --git a/tests/Fist/sample_psf/bromochloroamine.psf b/tests/Fist/sample_psf/bromochloroamine.psf new file mode 100644 index 0000000000..9ef67e23b7 --- /dev/null +++ b/tests/Fist/sample_psf/bromochloroamine.psf @@ -0,0 +1,28 @@ +PSF + + 1 !NTITLE +None + + 4 !NATOM + 1 NM00 1 NAME Br Br 0.0000 79.904000 0 + 2 NM00 1 NAME Cl Cl 0.0000 35.453000 0 + 3 NM00 1 NAME C C 0.0000 14.006700 0 + 4 NM00 1 NAME H H 0.0000 1.007940 0 + + 3 !NBOND + 3 1 3 2 4 3 + + 3 !NTHETA + 1 3 2 1 3 4 2 3 4 + + 0 !NPHI + + 1 !NIMPHI + 3 1 2 4 + + 0 !NDON + + 0 !NNB + + 0 !NGRP + diff --git a/tests/Fist/sample_psf/imp_test.psf b/tests/Fist/sample_psf/imp_test.psf index 3e0c619131..77169b371e 100644 --- a/tests/Fist/sample_psf/imp_test.psf +++ b/tests/Fist/sample_psf/imp_test.psf @@ -13,12 +13,12 @@ None 3 1 3 2 4 3 3 !NTHETA - 1 3 2 1 3 4 2 3 4 + 1 3 2 1 3 4 2 3 4 0 !NPHI - 1 !NIMPHI - 3 2 1 4 + 2 !NIMPHI + 3 1 2 4 3 2 1 4 0 !NDON