32 neighbor_list_iterator_p_type,&
34 neighbor_list_set_p_type
37 #include "./base/base_uses.f90"
43 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'qs_gcp_method'
64 TYPE(qs_environment_type),
POINTER :: qs_env
65 TYPE(qs_gcp_type),
POINTER :: gcp_env
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
89 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
90 TYPE(atprop_type),
POINTER :: atprop
91 TYPE(cell_type),
POINTER :: cell
92 TYPE(mp_para_env_type),
POINTER :: para_env
93 TYPE(neighbor_list_iterator_p_type), &
94 DIMENSION(:),
POINTER :: nl_iterator
95 TYPE(neighbor_list_set_p_type),
DIMENSION(:), &
97 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
98 TYPE(qs_force_type),
DIMENSION(:),
POINTER :: force
99 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
100 TYPE(virial_type),
POINTER :: virial
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)
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Calculation of the overlap integrals over Cartesian Gaussian-type functions.
subroutine, public overlap_ab(la_max, la_min, npgfa, rpgfa, zeta, lb_max, lb_min, npgfb, rpgfb, zetb, rab, sab, dab, ddab)
Calculation of the two-center overlap integrals [a|b] over Cartesian Gaussian-type functions....
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Holds information on atomic properties.
subroutine, public atprop_array_init(atarray, natom)
...
Handles all functions related to the CELL.
various routines to log and control the output. The idea is that decisions about where to log should ...
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Defines the basic variable types.
integer, parameter, public dp
Interface to the message passing library MPI.
Define the data structure for the particle information.
Definition of physical constants:
real(kind=dp), parameter, public kcalmol
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.
Calculation of gCP pair potentials.
subroutine, public calculate_gcp_pairpot(qs_env, gcp_env, energy, calculate_forces, ategcp)
...
Definition of gCP types for DFT calculations.
Define the quickstep kind type and their sub types.
Define the neighbor list data types and the corresponding functionality.
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
subroutine, public neighbor_list_iterator_release(iterator_set)
...
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
pure subroutine, public virial_pair_force(pv_virial, f0, force, rab)
Computes the contribution to the stress tensor from two-body pair-wise forces.