66 REAL(kind=
dp),
INTENT(OUT) :: energy
67 LOGICAL,
INTENT(IN) :: calculate_forces
68 REAL(kind=
dp),
DIMENSION(:),
OPTIONAL :: ategcp
70 CHARACTER(LEN=*),
PARAMETER :: routinen =
'calculate_gcp_pairpot'
72 INTEGER :: atom_a, atom_b, handle, i, iatom, ikind, &
73 jatom, jkind, mepos, natom, nkind, &
75 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: atom_of_kind, kind_of, ngcpat
76 LOGICAL :: atenergy, atex, atstress, use_virial, &
78 REAL(kind=
dp) :: eama, eamb, egcp, expab,
fac, fda, fdb, &
79 gnorm, nvirta, nvirtb, rcc, sint, sqa, &
81 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: egcpat
82 REAL(kind=
dp),
DIMENSION(3) :: dsint, fdij, rij
83 REAL(kind=
dp),
DIMENSION(3, 3) :: dvirial
84 REAL(kind=
dp),
DIMENSION(6) :: cla, clb, rcut, zeta, zetb
85 REAL(kind=
dp),
DIMENSION(6, 6) :: sab
86 REAL(kind=
dp),
DIMENSION(6, 6, 3) :: dab
87 REAL(kind=
dp),
DIMENSION(:),
POINTER :: atener
88 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: atstr
94 DIMENSION(:),
POINTER :: nl_iterator
103 IF (.NOT. gcp_env%do_gcp)
RETURN
105 CALL timeset(routinen, handle)
107 NULLIFY (atomic_kind_set, qs_kind_set, particle_set, sab_gcp)
109 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
110 cell=cell, virial=virial, para_env=para_env, atprop=atprop)
111 nkind =
SIZE(atomic_kind_set)
112 NULLIFY (particle_set)
113 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
114 natom =
SIZE(particle_set)
116 verbose = gcp_env%verbose
123 atenergy = atprop%energy
126 atener => atprop%ategcp
128 atstress = atprop%stress
129 atstr => atprop%atstress
132 IF (
PRESENT(ategcp)) atex = .true.
134 IF (unit_nr > 0)
THEN
136 WRITE (unit_nr, *)
" Pair potential geometrical counterpoise (gCP) calculation"
138 WRITE (unit_nr,
"(T15,A,T74,F7.4)")
" Gloabal Parameters: sigma = ", gcp_env%sigma, &
139 " alpha = ", gcp_env%alpha, &
140 " beta = ", gcp_env%beta, &
141 " eta = ", gcp_env%eta
143 WRITE (unit_nr,
"(T31,4(A5,10X))")
" kind",
"nvirt",
"Emiss",
" asto"
145 WRITE (unit_nr,
"(T31,i5,F15.1,F15.4,F15.4)") ikind, gcp_env%gcp_kind(ikind)%nbvirt, &
146 gcp_env%gcp_kind(ikind)%eamiss, gcp_env%gcp_kind(ikind)%asto
151 IF (calculate_forces)
THEN
154 ALLOCATE (atom_of_kind(natom), kind_of(natom))
156 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
157 IF (use_virial) dvirial = virial%pv_virial
165 ALLOCATE (egcpat(natom), ngcpat(natom))
172 cpassert(nsto ==
SIZE(gcp_env%gcp_kind(jkind)%al))
175 sab_gcp => gcp_env%sab_gcp
179 CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, r=rij)
181 rcc = sqrt(rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3))
182 IF (rcc > 1.e-6_dp)
THEN
184 IF (iatom == jatom)
fac = 0.5_dp
185 nvirta = gcp_env%gcp_kind(ikind)%nbvirt
186 nvirtb = gcp_env%gcp_kind(jkind)%nbvirt
187 eama = gcp_env%gcp_kind(ikind)%eamiss
188 eamb = gcp_env%gcp_kind(jkind)%eamiss
189 expab = exp(-gcp_env%alpha*rcc**gcp_env%beta)
190 zeta(1:nsto) = gcp_env%gcp_kind(ikind)%al(1:nsto)
191 zetb(1:nsto) = gcp_env%gcp_kind(jkind)%al(1:nsto)
192 cla(1:nsto) = gcp_env%gcp_kind(ikind)%cl(1:nsto)
193 clb(1:nsto) = gcp_env%gcp_kind(jkind)%cl(1:nsto)
194 IF (calculate_forces)
THEN
195 CALL overlap_ab(0, 0, nsto, rcut, zeta, 0, 0, nsto, rcut, zetb, rij, sab, dab)
197 dsint(i) = sum(cla*matmul(dab(:, :, i), clb))
200 CALL overlap_ab(0, 0, nsto, rcut, zeta, 0, 0, nsto, rcut, zetb, rij, sab)
202 sint = sum(cla*matmul(sab, clb))
203 IF (sint < 1.e-16_dp) cycle
204 sqa = sqrt(sint*nvirta)
205 sqb = sqrt(sint*nvirtb)
206 IF (sqb > 1.e-12_dp)
THEN
207 fda = gcp_env%sigma*eama*expab/sqb
211 IF (sqa > 1.e-12_dp)
THEN
212 fdb = gcp_env%sigma*eamb*expab/sqa
216 egcp = egcp +
fac*(fda + fdb)
218 egcpat(iatom) = egcpat(iatom) +
fac*fda
219 egcpat(jatom) = egcpat(jatom) +
fac*fdb
220 ngcpat(iatom) = ngcpat(iatom) + 1
221 ngcpat(jatom) = ngcpat(jatom) + 1
223 IF (calculate_forces)
THEN
224 fdij = -
fac*(fda + fdb)*(gcp_env%alpha*gcp_env%beta*rcc**(gcp_env%beta - 1.0_dp)*rij(1:3)/rcc)
225 IF (sqa > 1.e-12_dp)
THEN
226 fdij = fdij + 0.25_dp*
fac*fdb/(sqa*sqa)*dsint(1:3)
228 IF (sqb > 1.e-12_dp)
THEN
229 fdij = fdij + 0.25_dp*
fac*fda/(sqb*sqb)*dsint(1:3)
231 atom_a = atom_of_kind(iatom)
232 atom_b = atom_of_kind(jatom)
233 force(ikind)%gcp(:, atom_a) = force(ikind)%gcp(:, atom_a) - fdij(:)
234 force(jkind)%gcp(:, atom_b) = force(jkind)%gcp(:, atom_b) + fdij(:)
244 atener(iatom) = atener(iatom) + fda*
fac
245 atener(jatom) = atener(jatom) + fdb*
fac
248 ategcp(iatom) = ategcp(iatom) + fda*
fac
249 ategcp(jatom) = ategcp(jatom) + fdb*
fac
257 CALL para_env%sum(egcp)
260 CALL para_env%sum(egcpat)
261 CALL para_env%sum(ngcpat)
264 IF (unit_nr > 0)
THEN
265 WRITE (unit_nr,
"(T15,A,T61,F20.10)")
" Total gCP energy [au] :", egcp
266 WRITE (unit_nr,
"(T15,A,T61,F20.10)")
" Total gCP energy [kcal] :", egcp*
kcalmol
268 WRITE (unit_nr,
"(T19,A)")
" gCP atomic energy contributions"
269 WRITE (unit_nr,
"(T19,A,T60,A20)")
" # sites",
" BSSE [kcal/mol]"
271 WRITE (unit_nr,
"(12X,I8,10X,I8,T61,F20.10)") i, ngcpat(i), egcpat(i)*
kcalmol
274 IF (calculate_forces)
THEN
275 IF (unit_nr > 0)
THEN
276 WRITE (unit_nr, *)
" gCP Forces "
277 WRITE (unit_nr, *)
" Atom Kind Forces "
281 ikind = kind_of(iatom)
282 atom_a = atom_of_kind(iatom)
283 fdij(1:3) = force(ikind)%gcp(:, atom_a)
284 CALL para_env%sum(fdij)
285 gnorm = gnorm + sum(abs(fdij))
286 IF (unit_nr > 0)
WRITE (unit_nr,
"(i5,i7,3F20.14)") iatom, ikind, fdij
288 IF (unit_nr > 0)
THEN
290 WRITE (unit_nr, *)
" |G| = ", gnorm
294 dvirial = virial%pv_virial - dvirial
295 CALL para_env%sum(dvirial)
296 IF (unit_nr > 0)
THEN
297 WRITE (unit_nr, *)
" Stress Tensor (gCP)"
298 WRITE (unit_nr,
"(3G20.12)") dvirial
299 WRITE (unit_nr, *)
" Tr(P)/3 : ", (dvirial(1, 1) + dvirial(2, 2) + dvirial(3, 3))/3._dp
305 DEALLOCATE (egcpat, ngcpat)
308 CALL timestop(handle)
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.