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

out of plane bends + regtests (Louis)

Browse source 
svn-origin-rev: 9215
This commit is contained in:
Teodoro Laino 2009-11-23 14:22:14 +00:00
parent 42f66e7679
commit ed13711e52
16 changed files with 1377 additions and 277 deletions
Download patch file Download diff file Expand all files Collapse all files

76
src/fist_force.F
View file

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

143
src/fist_intra_force.F
View file

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

99
src/force_field_types.F
View file

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

315
src/force_fields_all.F
View file

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

60
src/force_fields_input.F
View file

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

177
src/force_fields_util.F
View file

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

170
src/input_constants.F
View file

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

75
src/input_cp2k_mm.F
View file

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

137
src/mol_force.F
View file

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

104
src/molecule_kind_types.F
View file

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

80
src/topology_connectivity_util.F
View file

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

1
tests/Fist/regtest-14/TEST_FILES
View file

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

5
tests/Fist/regtest-14/imp_test_12.inp
View file

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

178
tests/Fist/regtest-14/opbend_test.inp Normal file
View file

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

28
tests/Fist/sample_psf/bromochloroamine.psf Normal file
View file

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

6
tests/Fist/sample_psf/imp_test.psf
View file

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