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

gCP potentials and Short range bond correction for vdW D3 method (general formula) (#417)

Browse source 
* Short range bond correction for vdW D3 method, parameterize general formula

* gCP potentials input

* gCP energy/forces

* Debug forces

* gCP method

* Update to latest trunk version

* Adjust input interface

* Copyright banner and remove replicated file
This commit is contained in:
Juerg Hutter 2019-06-17 17:07:29 +02:00 • committed by GitHub
parent d22e34f677
commit ec01bc2fbd
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
29 changed files with 3315 additions and 20 deletions
Download patch file Download diff file Expand all files Collapse all files

1703
data/BASIS_MINIX Normal file
View file

File diff suppressed because it is too large Load diff

62
data/xc_section/HF-3c.sec Normal file
View file

@ -0,0 +1,62 @@
&XC
&XC_FUNCTIONAL NONE
&END XC_FUNCTIONAL
&HF
&SCREENING
EPS_SCHWARZ 1.0E-12
&END
&END
&vdW_POTENTIAL
DISPERSION_FUNCTIONAL PAIR_POTENTIAL
&PAIR_POTENTIAL
TYPE DFTD3(BJ)
D3BJ_SCALING 1.0000 0.4171 0.8777 2.9149
PARAMETER_FILE_NAME dftd3.dat
SHORT_RANGE_CORRECTION TRUE
SHORT_RANGE_CORRECTION_PARAMETERS 0.03 0.70 1.50 0.75
&END PAIR_POTENTIAL
&END vdW_POTENTIAL
&gcp_potential
GLOBAL_PARAMETERS 0.1290 1.1549 1.1763 1.1526
DELTA_ENERGY H 0.04240
DELTA_ENERGY He 0.02832
DELTA_ENERGY Li 0.17787
DELTA_ENERGY Be 0.17160
DELTA_ENERGY B 0.22424
DELTA_ENERGY C 0.27995
DELTA_ENERGY N 0.35791
DELTA_ENERGY O 0.47901
DELTA_ENERGY F 0.63852
DELTA_ENERGY Ne 0.83235
DELTA_ENERGY Na 1.11411
DELTA_ENERGY Mg 1.27115
DELTA_ENERGY Al 1.44695
DELTA_ENERGY Si 1.61098
DELTA_ENERGY P 1.76661
DELTA_ENERGY S 1.98823
DELTA_ENERGY Cl 2.22845
DELTA_ENERGY Ar 2.48796
DELTA_ENERGY K 0.37425
DELTA_ENERGY Ca 0.46097
DELTA_ENERGY Sc 0.44489
DELTA_ENERGY Ti 0.40499
DELTA_ENERGY V 0.37841
DELTA_ENERGY Cr 0.37344
DELTA_ENERGY Mn 0.36125
DELTA_ENERGY Fe 0.36001
DELTA_ENERGY Co 0.36293
DELTA_ENERGY Ni 0.24380
DELTA_ENERGY Cu 0.40530
DELTA_ENERGY Zn 0.39651
DELTA_ENERGY Ga 0.35002
DELTA_ENERGY Ge 0.34578
DELTA_ENERGY As 0.34953
DELTA_ENERGY Se 0.36731
DELTA_ENERGY Br 0.38201
DELTA_ENERGY Kr 0.39971
&END GCP_POTENTIAL
&END XC

19
src/common/bibliography.F
View file

@ -69,7 +69,7 @@ MODULE bibliography
Putrino2002, Sebastiani2001, Weber2009, Tran2013, Golze2013, &
Golze2015, Golze2017a, Golze2017b, Tuckerman1992, Zhao1994, &
Tozer1996, Goedecker2004, &
Khaliullin2013, Hutter2014, Bengtsson1999, Kantorovich2008, &
Khaliullin2013, Kruse2012, Hutter2014, Bengtsson1999, Kantorovich2008, &
Kantorovich2008a, Wellendorff2012, Niklasson2014, Borstnik2014, &
Rayson2009, Grimme2011, Fattebert2002, Andreussi2012, &
Khaliullin2007, Khaliullin2008, Merlot2014, Lin2009, Lin2013, &
@ -3983,6 +3983,23 @@ CONTAINS
"ER"), &
DOI="10.4208/cicp.OA-2018-0053")
CALL add_reference(key=Kruse2012, ISI_record=s2a( &
"AU Kruse,Holger", &
" Grimme,Stefan", &
"AF Kruse,Holger", &
" Grimme,Stefan", &
"TI A geometrical correction for the inter- and intra-molecular basis set superposition error", &
" in Hartree-Fock and density functional theory calculations for large systems", &
"SO The Journal of Chemical Physics", &
"SN 0021-9606", &
"PY 2012", &
"VL 136", &
"IS 15", &
"AR 154101", &
"DI 10.1063/1.3700154", &
"ER"), &
DOI="10.1063/1.3700154")
END SUBROUTINE add_all_references
END MODULE bibliography

68
src/input_cp2k_xc.F
View file

@ -12,7 +12,7 @@
MODULE input_cp2k_xc
USE bibliography, ONLY: &
Becke1988, Becke1997, BeckeRoussel1989, Goedecker1996, Grimme2006, Grimme2010, Grimme2011, &
Heyd2004, Lee1988, Lehtola2018, Marques2012, Ortiz1994, Perdew1981, Perdew1996, &
Heyd2004, Kruse2012, Lee1988, Lehtola2018, Marques2012, Ortiz1994, Perdew1981, Perdew1996, &
Perdew2008, Proynov2007, Tao2003, Tran2013, Vosko1980, Wellendorff2012, Zhang1998
USE cp_output_handling, ONLY: add_last_numeric,&
cp_print_key_section_create,&
@ -961,8 +961,16 @@ CONTAINS
usage="SHORT_RANGE_CORRECTION", default_l_val=.FALSE., &
lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(subsection, keyword)
! KG molecular corrections
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="SHORT_RANGE_CORRECTION_PARAMETERS", &
description="Parameters for the short-range bond correction to the DFT-D3 model. "// &
" s*(za*zb)^t1*EXP(-g*dr*r0ab^t2), parameters: s, g, t1, t2 "// &
" Defaults: s=0.08, g=10.0, t1=0.5, t2=-1.0 ", &
usage="SHORT_RANGE_CORRECTION_PARAMETRS", default_r_vals=(/0.08_dp, 10.0_dp, 0.5_dp, -1.0_dp/), &
n_var=4, type_of_var=real_t)
CALL section_add_keyword(subsection, keyword)
CALL keyword_release(keyword)
! KG molecular corrections
CALL keyword_create(keyword, __LOCATION__, name="MOLECULE_CORRECTION", &
description="Calculate a intermolecular correction to the DFT-D3 model", &
usage="MOLECULE_CORRECTION", default_l_val=.FALSE., &
@ -1075,6 +1083,58 @@ CONTAINS
END SUBROUTINE create_vdw_potential_section
! **************************************************************************************************
!> \brief creates the structure of the section needed for gCP potentials
!> \param section the section that will be created
!> \author jgh
! **************************************************************************************************
SUBROUTINE create_gcp_potential_section(section)
TYPE(section_type), POINTER :: section
CHARACTER(len=*), PARAMETER :: routineN = 'create_gcp_potential_section', &
routineP = moduleN//':'//routineN
TYPE(keyword_type), POINTER :: keyword
CPASSERT(.NOT. ASSOCIATED(section))
CALL section_create(section, __LOCATION__, name="gcp_potential", &
description="This section combines geometrical counterpoise potentials."// &
" This is a simple empirical pair potential to correct for BSSE. ", &
citations=(/Kruse2012/), &
n_keywords=1, n_subsections=1, repeats=.FALSE.)
NULLIFY (keyword)
CALL keyword_create(keyword, __LOCATION__, name="PARAMETER_FILE_NAME", &
description="Name of the parameter file, may include a path", &
usage="PARAMETER_FILE_NAME <FILENAME>", &
default_lc_val="---")
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="GLOBAL_PARAMETERS", &
description="Global parameters of the gCP method."// &
" Parameters are sigma, alpha, beta, eta from the original paper.", &
usage="GLOBAL_PARAMETERS 1.0 1.0 1.0 1.0", n_var=4, &
default_r_vals=(/0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp/))
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="DELTA_ENERGY", &
description="Specify the delta energy [Hartree] term for an atom kind", &
usage="DELTA_ENERGY type value", &
type_of_var=char_t, repeats=.TRUE., n_var=-1, default_c_vals=(/"XX ", "0.0"/))
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="VERBOSE", &
description="Verbose output for gCP calculation", &
usage="VERBOSE logical_value", &
default_l_val=.TRUE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
END SUBROUTINE create_gcp_potential_section
! **************************************************************************************************
!> \brief creates the input section for the xc part
!> \param section the section to create
@ -1196,6 +1256,10 @@ CONTAINS
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
CALL create_gcp_potential_section(subsection)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
END SUBROUTINE create_xc_section
! **************************************************************************************************

3
src/kg_energy_corrections.F
View file

@ -776,6 +776,9 @@ CONTAINS
! dispersion through pairpotentials
CALL calculate_dispersion_pairpot(qs_env, ec_env%dispersion_env, energy, .FALSE.)
ec_env%edispersion = ec_env%edispersion+energy
!deb ! gCP pairpotentials
!deb CALL calculate_gcp_pairpot(qs_env, ec_env%gcp_env, energy, .FALSE.)
!deb ec_env%edispersion = ec_env%edispersion+energy
SELECT CASE (ec_env%energy_functional)
CASE (kg_ec_functional_harris)

14
src/qs_dispersion_pairpot.F
View file

@ -873,8 +873,8 @@ CONTAINS
dc8a, dc8b, dcc6aba, dcc6abb, dcc6bcb, dcc6bcc, dcc6caa, dcc6cac, dd, de6, de8, de91, &
de921, de922, dea, devdw, dfdab6, dfdab8, dfdabc, dfdmp, dr, drk, e6, e6tot, e8, e8tot, &
e9, e9tot, elrc6, elrc8, elrc9, eps_cn, er, esrb, evdw, f0ab, fac, fac0, fdab6, fdab8, &
fdabc, fdmp, gnorm, kgc8, nab, nabc, r0, r0ab, r2ab, r2abc, r2bc, r2ca, r6, r8, rabc, &
rc2, rcc, rcut, rs6, rs8, s1, s2, s3, s6, s8, s8i, s9, srbe, xp
fdabc, fdmp, gnorm, gsrb, kgc8, nab, nabc, r0, r0ab, r2ab, r2abc, r2bc, r2ca, r6, r8, &
rabc, rc2, rcc, rcut, rs6, rs8, s1, s2, s3, s6, s8, s8i, s9, srbe, ssrb, t1srb, t2srb, xp
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: atom2mol, c6d2, cnkind, cnumbers, &
cnumfix, radd2
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: rcpbc
@ -1035,6 +1035,11 @@ CONTAINS
ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN
idmp = 2
END IF
! SRB parameters
ssrb = dispersion_env%srb_params(1)
gsrb = dispersion_env%srb_params(2)
t1srb = dispersion_env%srb_params(3)
t2srb = dispersion_env%srb_params(4)
IF (unit_nr > 0) THEN
WRITE (unit_nr, *) " Scaling parameter (s6) ", s6
@ -1049,6 +1054,7 @@ CONTAINS
WRITE (unit_nr, *) " Cutoff coordination numbers", eps_cn
IF (dispersion_env%lrc) THEN
WRITE (unit_nr, *) " Apply a long range correction"
WRITE (unit_nr, *) " SRB parameters (s,g,t1,t2) ", ssrb, gsrb, t1srb, t2srb
END IF
IF (dispersion_env%srb) THEN
WRITE (unit_nr, *) " Apply a short range bond correction"
@ -1239,7 +1245,7 @@ CONTAINS
CPABORT("Unknown DFT-D3 damping function:")
END IF
IF (dispersion_env%srb .AND. dr .LT. 30.0d0) THEN
srbe = 0.08_dp*SQRT(REAL((za*zb), KIND=dp))*EXP(-10.0_dp*dr/dispersion_env%r0ab(za, zb))
srbe = ssrb*(REAL((za*zb), KIND=dp))**t1srb*EXP(-gsrb*dr*dispersion_env%r0ab(za, zb)**t2srb)
esrb = esrb+srbe
evdw = evdw-srbe
ELSE
@ -1267,7 +1273,7 @@ CONTAINS
END IF
fdij(:) = (de6+de8)*rij(:)/dr*fac
IF (dispersion_env%srb .AND. dr .LT. 30.0d0) THEN
fdij(:) = fdij(:)+10._dp*srbe*rij(:)/(dr*dispersion_env%r0ab(za, zb))
fdij(:) = fdij(:)+srbe*gsrb*dispersion_env%r0ab(za, zb)**t2srb*rij(:)/dr
END IF
atom_a = atom_of_kind(iatom)
atom_b = atom_of_kind(jatom)

1
src/qs_dispersion_types.F
View file

@ -50,6 +50,7 @@ MODULE qs_dispersion_types
LOGICAL :: c9cnst !use constant c9 terms
LOGICAL :: lrc !calculate a long range correction
LOGICAL :: srb !calculate a short range bond correction
REAL(KIND=dp), DIMENSION(4) :: srb_params ! parameters for SRB (s,g,t1,t2)
TYPE(neighbor_list_set_p_type), &
DIMENSION(:), POINTER :: sab_vdw, sab_cn ! neighborlists for pair interactions
REAL(KIND=dp), DIMENSION(:, :, :, :, :), POINTER &

2
src/qs_dispersion_utils.F
View file

@ -113,6 +113,8 @@ CONTAINS
CALL section_vals_val_get(pp_section, "REFERENCE_C9_TERM", l_val=dispersion_env%c9cnst)
CALL section_vals_val_get(pp_section, "LONG_RANGE_CORRECTION", l_val=dispersion_env%lrc)
CALL section_vals_val_get(pp_section, "SHORT_RANGE_CORRECTION", l_val=dispersion_env%srb)
CALL section_vals_val_get(pp_section, "SHORT_RANGE_CORRECTION_PARAMETERS", r_vals=params)
dispersion_env%srb_params(1:4) = params(1:4)
! KG corrections
CALL section_vals_val_get(pp_section, "MOLECULE_CORRECTION", l_val=dispersion_env%domol)
CALL section_vals_val_get(pp_section, "MOLECULE_CORRECTION_C8", r_val=dispersion_env%kgc8)

9
src/qs_energy_init.F
View file

@ -48,6 +48,8 @@ MODULE qs_energy_init
USE qs_external_density, ONLY: external_read_density
USE qs_external_potential, ONLY: external_c_potential,&
external_e_potential
USE qs_gcp_method, ONLY: calculate_gcp_pairpot
USE qs_gcp_types, ONLY: qs_gcp_type
USE qs_ks_methods, ONLY: qs_ks_allocate_basics
USE qs_ks_types, ONLY: qs_ks_env_type,&
set_ks_env
@ -261,6 +263,7 @@ CONTAINS
POINTER :: sab_nl
TYPE(qs_dispersion_type), POINTER :: dispersion_env
TYPE(qs_energy_type), POINTER :: energy
TYPE(qs_gcp_type), POINTER :: gcp_env
TYPE(section_vals_type), POINTER :: input
CALL timeset(routineN, handle)
@ -339,6 +342,12 @@ CONTAINS
! energy info at the end of the SCF
CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env, energy=energy)
CALL calculate_dispersion_pairpot(qs_env, dispersion_env, energy%dispersion, calc_forces)
! Add possible pair potential gCP energy - Evaluate first so we can print
! energy info at the end of the SCF
CALL get_qs_env(qs_env=qs_env, gcp_env=gcp_env, energy=energy)
IF (ASSOCIATED(gcp_env)) THEN
CALL calculate_gcp_pairpot(qs_env, gcp_env, energy%gcp, calc_forces)
END IF
END IF
! Embedding potential

2
src/qs_energy_types.F
View file

@ -27,6 +27,7 @@ MODULE qs_energy_types
core_self, &
repulsive, &
dispersion, &
gcp, &
ex, &
exc, &
exc_aux_fit, &
@ -154,6 +155,7 @@ CONTAINS
qs_energy%core_self = 0.0_dp
qs_energy%repulsive = 0.0_dp
qs_energy%dispersion = 0.0_dp
qs_energy%gcp = 0.0_dp
qs_energy%qmmm_el = 0.0_dp
qs_energy%qmmm_nu = 0.0_dp
qs_energy%ex = 0.0_dp

1
src/qs_energy_utils.F
View file

@ -250,6 +250,7 @@ CONTAINS
CALL atprop_array_add(atprop%atener, atprop%atexc)
CALL atprop_array_add(atprop%atener, atprop%atecoul)
CALL atprop_array_add(atprop%atener, atprop%atevdw)
CALL atprop_array_add(atprop%atener, atprop%ategcp)
CALL atprop_array_add(atprop%atener, atprop%atecc)
CALL atprop_array_add(atprop%atener, atprop%ate1c)
END IF

19
src/qs_environment.F
View file

@ -142,6 +142,9 @@ MODULE qs_environment
qs_environment_type,&
set_qs_env
USE qs_force_types, ONLY: qs_force_type
USE qs_gcp_types, ONLY: qs_gcp_type
USE qs_gcp_utils, ONLY: qs_gcp_env_set,&
qs_gcp_init
USE qs_interactions, ONLY: init_interaction_radii,&
init_se_nlradius,&
write_core_charge_radii,&
@ -516,6 +519,7 @@ CONTAINS
TYPE(qs_dispersion_type), POINTER :: dispersion_env
TYPE(qs_energy_type), POINTER :: energy
TYPE(qs_force_type), DIMENSION(:), POINTER :: force
TYPE(qs_gcp_type), POINTER :: gcp_env
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_kind_type), POINTER :: qs_kind
TYPE(qs_ks_env_type), POINTER :: ks_env
@ -1141,6 +1145,21 @@ CONTAINS
CALL set_qs_env(qs_env, dispersion_env=dispersion_env)
END IF
! Initialize possible geomertical counterpoise correction potential
IF (dft_control%qs_control%method_id == do_method_gpw .OR. &
dft_control%qs_control%method_id == do_method_gapw .OR. &
dft_control%qs_control%method_id == do_method_gapw_xc .OR. &
dft_control%qs_control%method_id == do_method_lrigpw .OR. &
dft_control%qs_control%method_id == do_method_rigpw .OR. &
dft_control%qs_control%method_id == do_method_ofgpw) THEN
ALLOCATE (gcp_env)
NULLIFY (xc_section)
xc_section => section_vals_get_subs_vals(dft_section, "XC")
CALL qs_gcp_env_set(gcp_env, xc_section)
CALL qs_gcp_init(qs_env, gcp_env)
CALL set_qs_env(qs_env, gcp_env=gcp_env)
END IF
! *** Allocate the MO data types ***
CALL get_qs_kind_set(qs_kind_set, nsgf=n_ao, nelectron=nelectron)

20
src/qs_environment_types.F
View file

@ -94,6 +94,8 @@ MODULE qs_environment_types
qs_dispersion_type
USE qs_energy_types, ONLY: qs_energy_type
USE qs_force_types, ONLY: qs_force_type
USE qs_gcp_types, ONLY: qs_gcp_release,&
qs_gcp_type
USE qs_kind_types, ONLY: qs_kind_type
USE qs_ks_qmmm_types, ONLY: qs_ks_qmmm_env_type,&
qs_ks_qmmm_release,&
@ -277,6 +279,8 @@ MODULE qs_environment_types
TYPE(lri_density_type), POINTER :: lri_density
! Empirical dispersion
TYPE(qs_dispersion_type), POINTER :: dispersion_env
! Empirical geometrical BSSE correction
TYPE(qs_gcp_type), POINTER :: gcp_env
! Semi-empirical and DFTB types
TYPE(ewald_environment_type), POINTER :: ewald_env
TYPE(ewald_pw_type), POINTER :: ewald_pw
@ -456,6 +460,7 @@ CONTAINS
!> \param lri_env ...
!> \param lri_density ...
!> \param dispersion_env ...
!> \param gcp_env ...
!> \param vee ...
!> \param rho_external ...
!> \param external_vxc ...
@ -507,7 +512,7 @@ CONTAINS
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, admm_dm, &
lri_env, lri_density, dispersion_env, vee, rho_external, external_vxc, mask, &
lri_env, lri_density, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, &
mp2_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, rhs)
@ -621,6 +626,7 @@ CONTAINS
TYPE(lri_environment_type), OPTIONAL, POINTER :: lri_env
TYPE(lri_density_type), OPTIONAL, POINTER :: lri_density
TYPE(qs_dispersion_type), OPTIONAL, POINTER :: dispersion_env
TYPE(qs_gcp_type), OPTIONAL, POINTER :: gcp_env
TYPE(pw_p_type), OPTIONAL, POINTER :: vee
TYPE(qs_rho_type), OPTIONAL, POINTER :: rho_external
TYPE(pw_p_type), OPTIONAL, POINTER :: external_vxc, mask
@ -694,6 +700,7 @@ CONTAINS
IF (PRESENT(lri_env)) lri_env => qs_env%lri_env
IF (PRESENT(lri_density)) lri_density => qs_env%lri_density
IF (PRESENT(dispersion_env)) dispersion_env => qs_env%dispersion_env
IF (PRESENT(gcp_env)) gcp_env => qs_env%gcp_env
IF (PRESENT(run_rtp)) run_rtp = qs_env%run_rtp
IF (PRESENT(rtp)) rtp => qs_env%rtp
IF (PRESENT(ls_scf_env)) ls_scf_env => qs_env%ls_scf_env
@ -909,7 +916,7 @@ CONTAINS
NULLIFY (qs_env%efield)
NULLIFY (qs_env%lri_env)
NULLIFY (qs_env%lri_density)
NULLIFY (qs_env%dispersion_env)
NULLIFY (qs_env%gcp_env)
NULLIFY (qs_env%rtp)
NULLIFY (qs_env%mp2_env)
NULLIFY (qs_env%kg_env)
@ -1012,6 +1019,7 @@ CONTAINS
!> \param lri_env ...
!> \param lri_density ...
!> \param dispersion_env ...
!> \param gcp_env ...
!> \param mp2_env ...
!> \param kg_env ...
!> \param kpoints ...
@ -1039,7 +1047,7 @@ CONTAINS
linres_control, xas_env, cp_ddapc_env, cp_ddapc_ewald, &
outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, &
se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, ls_scf_env, &
do_transport, transport_env, lri_env, lri_density, dispersion_env, mp2_env, kg_env, &
do_transport, transport_env, lri_env, lri_density, dispersion_env, gcp_env, mp2_env, kg_env, &
kpoints, WannierCentres, almo_scf_env, gradient_history, variable_history, embed_pot, &
spin_embed_pot, polar_env, rhs)
@ -1101,6 +1109,7 @@ CONTAINS
TYPE(lri_environment_type), OPTIONAL, POINTER :: lri_env
TYPE(lri_density_type), OPTIONAL, POINTER :: lri_density
TYPE(qs_dispersion_type), OPTIONAL, POINTER :: dispersion_env
TYPE(qs_gcp_type), OPTIONAL, POINTER :: gcp_env
TYPE(mp2_type), OPTIONAL, POINTER :: mp2_env
TYPE(kg_environment_type), OPTIONAL, POINTER :: kg_env
TYPE(kpoint_type), OPTIONAL, POINTER :: kpoints
@ -1282,6 +1291,7 @@ CONTAINS
IF (PRESENT(lri_env)) qs_env%lri_env => lri_env
IF (PRESENT(lri_density)) qs_env%lri_density => lri_density
IF (PRESENT(dispersion_env)) qs_env%dispersion_env => dispersion_env
IF (PRESENT(gcp_env)) qs_env%gcp_env => gcp_env
IF (PRESENT(WannierCentres)) qs_env%WannierCentres => WannierCentres
IF (PRESENT(kpoints)) CALL set_ks_env(qs_env%ks_env, kpoints=kpoints)
@ -1510,6 +1520,10 @@ CONTAINS
! dispersion
CALL qs_dispersion_release(qs_env%dispersion_env)
! gCP
IF (ASSOCIATED(qs_env%gcp_env)) THEN
CALL qs_gcp_release(qs_env%gcp_env)
END IF
IF (ASSOCIATED(qs_env%WannierCentres)) THEN
DO i = 1, SIZE(qs_env%WannierCentres)

9
src/qs_force.F
View file

@ -552,9 +552,10 @@ CONTAINS
iatom, ikind, " core_overlap", qs_force(ikind)%core_overlap(1:3, i), &
iatom, ikind, " rho_core", qs_force(ikind)%rho_core(1:3, i), &
iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
iatom, ikind, " rho_lri_elec", qs_force(ikind)%rho_lri_elec(1:3, i), &
iatom, ikind, " rho_lri_elec", qs_force(ikind)%rho_lri_elec(1:3, i), &
iatom, ikind, " ch_pulay", qs_force(ikind)%ch_pulay(1:3, i), &
iatom, ikind, " dispersion", qs_force(ikind)%dispersion(1:3, i), &
iatom, ikind, " gCP", qs_force(ikind)%gcp(1:3, i), &
iatom, ikind, " other", qs_force(ikind)%other(1:3, i), &
iatom, ikind, " fock_4c", qs_force(ikind)%fock_4c(1:3, i), &
iatom, ikind, " ehrenfest", qs_force(ikind)%ehrenfest(1:3, i), &
@ -580,11 +581,12 @@ CONTAINS
iatom, ikind, " core_overlap", qs_force(ikind)%core_overlap(1:3, i), &
iatom, ikind, " rho_core", qs_force(ikind)%rho_core(1:3, i), &
iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
iatom, ikind, " rho_lri_elec", qs_force(ikind)%rho_lri_elec(1:3, i), &
iatom, ikind, " rho_lri_elec", qs_force(ikind)%rho_lri_elec(1:3, i), &
iatom, ikind, " vhxc_atom", qs_force(ikind)%vhxc_atom(1:3, i), &
iatom, ikind, " g0s_Vh_elec", qs_force(ikind)%g0s_Vh_elec(1:3, i), &
iatom, ikind, " ch_pulay", qs_force(ikind)%ch_pulay(1:3, i), &
iatom, ikind, " dispersion", qs_force(ikind)%dispersion(1:3, i), &
iatom, ikind, " gCP", qs_force(ikind)%gcp(1:3, i), &
iatom, ikind, " fock_4c", qs_force(ikind)%fock_4c(1:3, i), &
iatom, ikind, " efield", qs_force(ikind)%efield(1:3, i), &
iatom, ikind, " eev", qs_force(ikind)%eev(1:3, i), &
@ -600,7 +602,6 @@ CONTAINS
WRITE (UNIT=output_unit, FMT=fmtstr2) &
iatom, ikind, " all_potential", qs_force(ikind)%all_potential(1:3, i), &
iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
iatom, ikind, " rho_lri_elec", qs_force(ikind)%rho_lri_elec(1:3, i), &
iatom, ikind, " total", qs_force(ikind)%total(1:3, i)
grand_total(1:3) = grand_total(1:3)+qs_force(ikind)%total(1:3, i)
END DO
@ -618,7 +619,6 @@ CONTAINS
iatom, ikind, " core_overlap", qs_force(ikind)%core_overlap(1:3, i), &
iatom, ikind, " rho_core", qs_force(ikind)%rho_core(1:3, i), &
iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
iatom, ikind, " rho_lri_elec", qs_force(ikind)%rho_lri_elec(1:3, i), &
iatom, ikind, " ch_pulay", qs_force(ikind)%ch_pulay(1:3, i), &
iatom, ikind, " fock_4c", qs_force(ikind)%fock_4c(1:3, i), &
iatom, ikind, " mp2_non_sep", qs_force(ikind)%mp2_non_sep(1:3, i), &
@ -649,7 +649,6 @@ CONTAINS
iatom, ikind, " overlap", qs_force(ikind)%overlap(1:3, i), &
iatom, ikind, " kinetic", qs_force(ikind)%kinetic(1:3, i), &
iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
iatom, ikind, " rho_lri_elec", qs_force(ikind)%rho_lri_elec(1:3, i), &
iatom, ikind, " dispersion", qs_force(ikind)%dispersion(1:3, i), &
iatom, ikind, " all potential", qs_force(ikind)%all_potential(1:3, i), &
iatom, ikind, " other", qs_force(ikind)%other(1:3, i), &

9
src/qs_force_types.F
View file

@ -38,6 +38,7 @@ MODULE qs_force_types
g0s_Vh_elec, &
repulsive, &
dispersion, &
gcp, &
other, &
ch_pulay, &
fock_4c, &
@ -100,6 +101,7 @@ CONTAINS
ALLOCATE (qs_force(ikind)%g0s_Vh_elec(3, n))
ALLOCATE (qs_force(ikind)%repulsive(3, n))
ALLOCATE (qs_force(ikind)%dispersion(3, n))
ALLOCATE (qs_force(ikind)%gcp(3, n))
ALLOCATE (qs_force(ikind)%other(3, n))
ALLOCATE (qs_force(ikind)%ch_pulay(3, n))
ALLOCATE (qs_force(ikind)%ehrenfest(3, n))
@ -196,6 +198,10 @@ CONTAINS
DEALLOCATE (qs_force(ikind)%dispersion)
END IF
IF (ASSOCIATED(qs_force(ikind)%gcp)) THEN
DEALLOCATE (qs_force(ikind)%gcp)
END IF
IF (ASSOCIATED(qs_force(ikind)%other)) THEN
DEALLOCATE (qs_force(ikind)%other)
END IF
@ -270,6 +276,7 @@ CONTAINS
qs_force(ikind)%g0s_Vh_elec(:, :) = 0.0_dp
qs_force(ikind)%repulsive(:, :) = 0.0_dp
qs_force(ikind)%dispersion(:, :) = 0.0_dp
qs_force(ikind)%gcp(:, :) = 0.0_dp
qs_force(ikind)%other(:, :) = 0.0_dp
qs_force(ikind)%fock_4c(:, :) = 0.0_dp
qs_force(ikind)%ehrenfest(:, :) = 0.0_dp
@ -320,6 +327,7 @@ CONTAINS
CALL mp_sum(qs_force(ikind)%mp2_sep, para_env%group)
CALL mp_sum(qs_force(ikind)%repulsive, para_env%group)
CALL mp_sum(qs_force(ikind)%dispersion, para_env%group)
CALL mp_sum(qs_force(ikind)%gcp, para_env%group)
CALL mp_sum(qs_force(ikind)%ehrenfest, para_env%group)
qs_force(ikind)%total(:, :) = qs_force(ikind)%total(:, :)+ &
@ -341,6 +349,7 @@ CONTAINS
qs_force(ikind)%mp2_sep(:, :)+ &
qs_force(ikind)%repulsive(:, :)+ &
qs_force(ikind)%dispersion(:, :)+ &
qs_force(ikind)%gcp(:, :)+ &
qs_force(ikind)%ehrenfest(:, :)+ &
qs_force(ikind)%efield(:, :)+ &
qs_force(ikind)%eev(:, :)

306
src/qs_gcp_method.F Normal file
View file

@ -0,0 +1,306 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 - 2019 CP2K developers group !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \brief Calculation of gCP pair potentials
!> \author JGH
! **************************************************************************************************
MODULE qs_gcp_method
USE ai_overlap, ONLY: overlap_ab
USE atomic_kind_types, ONLY: atomic_kind_type,&
get_atomic_kind_set
USE atprop_types, ONLY: atprop_array_init,&
atprop_type
USE cell_types, ONLY: cell_type
USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
USE cp_para_types, ONLY: cp_para_env_type
USE kinds, ONLY: dp
USE message_passing, ONLY: mp_sum
USE particle_types, ONLY: particle_type
USE physcon, ONLY: kcalmol
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_force_types, ONLY: qs_force_type
USE qs_gcp_types, ONLY: qs_gcp_type
USE qs_kind_types, ONLY: qs_kind_type
USE qs_neighbor_list_types, ONLY: get_iterator_info,&
neighbor_list_iterate,&
neighbor_list_iterator_create,&
neighbor_list_iterator_p_type,&
neighbor_list_iterator_release,&
neighbor_list_set_p_type
USE virial_methods, ONLY: virial_pair_force
USE virial_types, ONLY: virial_type
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_gcp_method'
PUBLIC :: calculate_gcp_pairpot
! **************************************************************************************************
CONTAINS
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param gcp_env ...
!> \param energy ...
!> \param calculate_forces ...
!> \note
!> \note energy_correction_type: also add gcp_env and egcp to the type
!> \note
! **************************************************************************************************
SUBROUTINE calculate_gcp_pairpot(qs_env, gcp_env, energy, calculate_forces)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(qs_gcp_type), POINTER :: gcp_env
REAL(KIND=dp), INTENT(OUT) :: energy
LOGICAL, INTENT(IN) :: calculate_forces
CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_gcp_pairpot', &
routineP = moduleN//':'//routineN
INTEGER :: atom_a, atom_b, handle, i, iatom, ikind, &
jatom, jkind, mepos, natom, nkind, &
nsto, unit_nr
INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of, ngcpat
LOGICAL :: atenergy, atstress, use_virial, verbose
REAL(KIND=dp) :: eama, eamb, egcp, expab, fac, fda, fdb, &
gnorm, nvirta, nvirtb, rcc, sint, sqa, &
sqb
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: egcpat
REAL(KIND=dp), DIMENSION(3) :: dsint, fdij, rij
REAL(KIND=dp), DIMENSION(3, 3) :: dvirial
REAL(KIND=dp), DIMENSION(6) :: cla, clb, rcut, zeta, zetb
REAL(KIND=dp), DIMENSION(6, 6) :: sab
REAL(KIND=dp), DIMENSION(6, 6, 3) :: dab
REAL(KIND=dp), DIMENSION(:), POINTER :: atener
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: atstr
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(atprop_type), POINTER :: atprop
TYPE(cell_type), POINTER :: cell
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(neighbor_list_iterator_p_type), &
DIMENSION(:), POINTER :: nl_iterator
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_gcp
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_force_type), DIMENSION(:), POINTER :: force
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(virial_type), POINTER :: virial
energy = 0._dp
IF (.NOT. gcp_env%do_gcp) RETURN
CALL timeset(routineN, handle)
NULLIFY (atomic_kind_set, qs_kind_set, particle_set, sab_gcp)
CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
cell=cell, virial=virial, para_env=para_env, atprop=atprop)
nkind = SIZE(atomic_kind_set)
NULLIFY (particle_set)
CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
natom = SIZE(particle_set)
verbose = gcp_env%verbose
IF (verbose) THEN
unit_nr = cp_logger_get_default_io_unit()
ELSE
unit_nr = -1
END IF
! atomic energy and stress arrays
atenergy = atprop%energy
IF (atenergy) THEN
CALL atprop_array_init(atprop%ategcp, natom)
atener => atprop%ategcp
END IF
atstress = atprop%stress
IF (atstress) THEN
atstr => atprop%atstress
END IF
IF (unit_nr > 0) THEN
WRITE (unit_nr, *)
WRITE (unit_nr, *) " Pair potential geometrical counterpoise (gCP) calculation"
WRITE (unit_nr, *)
WRITE (unit_nr, "(T15,A,T74,F7.4)") " Gloabal Parameters: sigma = ", gcp_env%sigma, &
" alpha = ", gcp_env%alpha, &
" beta = ", gcp_env%beta, &
" eta = ", gcp_env%eta
WRITE (unit_nr, *)
WRITE (unit_nr, "(T31,4(A5,10X))") " kind", "nvirt", "Emiss", " asto"
DO ikind = 1, nkind
WRITE (unit_nr, "(T31,i5,F15.1,F15.4,F15.4)") ikind, gcp_env%gcp_kind(ikind)%nbvirt, &
gcp_env%gcp_kind(ikind)%eamiss, gcp_env%gcp_kind(ikind)%asto
END DO
WRITE (unit_nr, *)
END IF
IF (calculate_forces) THEN
NULLIFY (force)
CALL get_qs_env(qs_env=qs_env, force=force)
ALLOCATE (atom_of_kind(natom), kind_of(natom))
CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
IF (use_virial) dvirial = virial%pv_virial
END IF
! include all integrals in the list
rcut = 1.e6_dp
egcp = 0.0_dp
IF (verbose) THEN
ALLOCATE (egcpat(natom), ngcpat(natom))
egcpat = 0.0_dp
ngcpat = 0
END IF
nsto = 6
DO ikind = 1, nkind
CPASSERT(nsto == SIZE(gcp_env%gcp_kind(jkind)%al))
END DO
sab_gcp => gcp_env%sab_gcp
CALL neighbor_list_iterator_create(nl_iterator, sab_gcp)
DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, r=rij)
rcc = SQRT(rij(1)*rij(1)+rij(2)*rij(2)+rij(3)*rij(3))
IF (rcc > 1.e-6_dp) THEN
fac = 1._dp
IF (iatom == jatom) fac = 0.5_dp
nvirta = gcp_env%gcp_kind(ikind)%nbvirt
nvirtb = gcp_env%gcp_kind(jkind)%nbvirt
eama = gcp_env%gcp_kind(ikind)%eamiss
eamb = gcp_env%gcp_kind(jkind)%eamiss
expab = EXP(-gcp_env%alpha*rcc**gcp_env%beta)
zeta(1:nsto) = gcp_env%gcp_kind(ikind)%al(1:nsto)
zetb(1:nsto) = gcp_env%gcp_kind(jkind)%al(1:nsto)
cla(1:nsto) = gcp_env%gcp_kind(ikind)%cl(1:nsto)
clb(1:nsto) = gcp_env%gcp_kind(jkind)%cl(1:nsto)
IF (calculate_forces) THEN
CALL overlap_ab(0, 0, nsto, rcut, zeta, 0, 0, nsto, rcut, zetb, rij, sab, dab)
DO i = 1, 3
dsint(i) = SUM(cla*MATMUL(dab(:, :, i), clb))
END DO
ELSE
CALL overlap_ab(0, 0, nsto, rcut, zeta, 0, 0, nsto, rcut, zetb, rij, sab)
END IF
sint = SUM(cla*MATMUL(sab, clb))
IF (sint < 1.e-16_dp) CYCLE
sqa = SQRT(sint*nvirta)
sqb = SQRT(sint*nvirtb)
IF (sqb > 1.e-12_dp) THEN
fda = gcp_env%sigma*eama*expab/sqb
ELSE
fda = 0.0_dp
END IF
IF (sqa > 1.e-12_dp) THEN
fdb = gcp_env%sigma*eamb*expab/sqa
ELSE
fdb = 0.0_dp
END IF
egcp = egcp+fac*(fda+fdb)
IF (verbose) THEN
egcpat(iatom) = egcpat(iatom)+fac*fda
egcpat(jatom) = egcpat(jatom)+fac*fdb
ngcpat(iatom) = ngcpat(iatom)+1
ngcpat(jatom) = ngcpat(jatom)+1
END IF
IF (calculate_forces) THEN
fdij = -fac*(fda+fdb)*(gcp_env%alpha*gcp_env%beta*rcc**(gcp_env%beta-1.0_dp)*rij(1:3)/rcc)
IF (sqa > 1.e-12_dp) THEN
fdij = fdij+0.25_dp*fac*fdb/(sqa*sqa)*dsint(1:3)
END IF
IF (sqb > 1.e-12_dp) THEN
fdij = fdij+0.25_dp*fac*fda/(sqb*sqb)*dsint(1:3)
END IF
atom_a = atom_of_kind(iatom)
atom_b = atom_of_kind(jatom)
force(ikind)%gcp(:, atom_a) = force(ikind)%gcp(:, atom_a)-fdij(:)
force(jkind)%gcp(:, atom_b) = force(jkind)%gcp(:, atom_b)+fdij(:)
IF (use_virial) THEN
CALL virial_pair_force(virial%pv_virial, -1._dp, fdij, rij)
END IF
IF (atstress) THEN
CALL virial_pair_force(atstr(:, :, iatom), -0.5_dp, fdij, rij)
CALL virial_pair_force(atstr(:, :, jatom), -0.5_dp, fdij, rij)
END IF
END IF
IF (atenergy) THEN
atener(iatom) = atener(iatom)+fda*fac
atener(jatom) = atener(jatom)+fdb*fac
END IF
END IF
END DO
CALL neighbor_list_iterator_release(nl_iterator)
! set gCP energy
CALL mp_sum(egcp, para_env%group)
energy = egcp
IF (verbose) THEN
CALL mp_sum(egcpat, para_env%group)
CALL mp_sum(ngcpat, para_env%group)
END IF
IF (unit_nr > 0) THEN
WRITE (unit_nr, "(T15,A,T61,F20.10)") " Total gCP energy [au] :", egcp
WRITE (unit_nr, "(T15,A,T61,F20.10)") " Total gCP energy [kcal] :", egcp*kcalmol
WRITE (unit_nr, *)
WRITE (unit_nr, "(T19,A)") " gCP atomic energy contributions"
WRITE (unit_nr, "(T19,A,T60,A20)") " # sites", " BSSE [kcal/mol]"
DO i = 1, natom
WRITE (unit_nr, "(12X,I8,10X,I8,T61,F20.10)") i, ngcpat(i), egcpat(i)*kcalmol
END DO
END IF
IF (calculate_forces) THEN
IF (unit_nr > 0) THEN
WRITE (unit_nr, *) " gCP Forces "
WRITE (unit_nr, *) " Atom Kind Forces "
END IF
gnorm = 0._dp
DO iatom = 1, natom
ikind = kind_of(iatom)
atom_a = atom_of_kind(iatom)
fdij(1:3) = force(ikind)%gcp(:, atom_a)
CALL mp_sum(fdij, para_env%group)
gnorm = gnorm+SUM(ABS(fdij))
IF (unit_nr > 0) WRITE (unit_nr, "(i5,i7,3F20.14)") iatom, ikind, fdij
END DO
IF (unit_nr > 0) THEN
WRITE (unit_nr, *)
WRITE (unit_nr, *) " |G| = ", gnorm
WRITE (unit_nr, *)
END IF
IF (use_virial) THEN
dvirial = virial%pv_virial-dvirial
CALL mp_sum(dvirial, para_env%group)
IF (unit_nr > 0) THEN
WRITE (unit_nr, *) " Stress Tensor (gCP)"
WRITE (unit_nr, "(3G20.12)") dvirial
WRITE (unit_nr, *) " Tr(P)/3 : ", (dvirial(1, 1)+dvirial(2, 2)+dvirial(3, 3))/3._dp
WRITE (unit_nr, *)
END IF
END IF
END IF
IF (verbose) THEN
DEALLOCATE (egcpat, ngcpat)
END IF
CALL timestop(handle)
END SUBROUTINE calculate_gcp_pairpot
! **************************************************************************************************
END MODULE qs_gcp_method

100
src/qs_gcp_types.F Normal file
View file

@ -0,0 +1,100 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 - 2019 CP2K developers group !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \brief Definition of gCP types for DFT calculations
!> \author JGH (20.10.2018)
! **************************************************************************************************
MODULE qs_gcp_types
USE kinds, ONLY: default_string_length,&
dp
USE qs_neighbor_list_types, ONLY: deallocate_neighbor_list_set,&
neighbor_list_set_p_type
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_gcp_types'
! **************************************************************************************************
TYPE qs_gcp_kind_type
INTEGER :: za
REAL(KIND=dp) :: asto
REAL(KIND=dp) :: rcsto
INTEGER :: nq
REAL(KIND=dp) :: nbvirt
REAL(KIND=dp) :: eamiss
REAL(KIND=dp), DIMENSION(6) :: al
REAL(KIND=dp), DIMENSION(6) :: cl
END TYPE qs_gcp_kind_type
! **************************************************************************************************
TYPE qs_gcp_type
LOGICAL :: do_gcp = .FALSE.
LOGICAL :: verbose !extended output
CHARACTER(LEN=default_string_length) :: parameter_file_name
! parameter input
CHARACTER(LEN=default_string_length), &
DIMENSION(:), POINTER :: kind_type => NULL()
REAL(KIND=dp), DIMENSION(:), POINTER :: ea => NULL()
!global parameters
REAL(KIND=dp) :: alpha, beta, sigma, eta
!neighborlist
TYPE(neighbor_list_set_p_type), &
DIMENSION(:), POINTER :: sab_gcp => Null() ! neighborlists for pair interactions
!kind information
TYPE(qs_gcp_kind_type), DIMENSION(:), &
POINTER :: gcp_kind => NULL() ! atomic kind parameters
END TYPE qs_gcp_type
! **************************************************************************************************
PUBLIC :: qs_gcp_type
PUBLIC :: qs_gcp_release
! **************************************************************************************************
CONTAINS
! **************************************************************************************************
!> \brief ...
!> \param gcp_env ...
! **************************************************************************************************
SUBROUTINE qs_gcp_release(gcp_env)
TYPE(qs_gcp_type), POINTER :: gcp_env
CHARACTER(len=*), PARAMETER :: routineN = 'qs_gcp_release', routineP = moduleN//':'//routineN
INTEGER :: iab
IF (ASSOCIATED(gcp_env)) THEN
IF (ASSOCIATED(gcp_env%sab_gcp)) THEN
DO iab = 1, SIZE(gcp_env%sab_gcp)
CALL deallocate_neighbor_list_set(gcp_env%sab_gcp(iab)%neighbor_list_set)
END DO
DEALLOCATE (gcp_env%sab_gcp)
END IF
IF (ASSOCIATED(gcp_env%kind_type)) THEN
DEALLOCATE (gcp_env%kind_type)
END IF
IF (ASSOCIATED(gcp_env%ea)) THEN
DEALLOCATE (gcp_env%ea)
END IF
IF (ASSOCIATED(gcp_env%gcp_kind)) THEN
DEALLOCATE (gcp_env%gcp_kind)
END IF
DEALLOCATE (gcp_env)
END IF
END SUBROUTINE qs_gcp_release
! **************************************************************************************************
END MODULE qs_gcp_types

251
src/qs_gcp_utils.F Normal file
View file

@ -0,0 +1,251 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 - 2019 CP2K developers group !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \brief Set disperson types for DFT calculations
!> \author JGH (04.2014)
! **************************************************************************************************
MODULE qs_gcp_utils
USE basis_set_types, ONLY: get_gto_basis_set,&
gto_basis_set_type
USE cp_para_types, ONLY: cp_para_env_type
USE cp_parser_methods, ONLY: parser_get_next_line,&
parser_get_object
USE cp_parser_types, ONLY: cp_parser_type,&
parser_create,&
parser_release
USE input_section_types, ONLY: section_vals_get,&
section_vals_get_subs_vals,&
section_vals_type,&
section_vals_val_get
USE kinds, ONLY: default_string_length,&
dp
USE mathconstants, ONLY: pi
USE periodic_table, ONLY: get_ptable_info
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_gcp_types, ONLY: qs_gcp_type
USE qs_kind_types, ONLY: get_qs_kind,&
qs_kind_type
USE sto_ng, ONLY: get_sto_ng
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_gcp_utils'
INTEGER, DIMENSION(106) :: nshell = (/ &
1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, & ! 1-30
4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, & ! 31-60
6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, & ! 61-90
7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7/) ! 91-106
INTEGER, DIMENSION(106) :: nll = (/ &
1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, & ! 1-30
2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 4, 4, 4, 4, & ! 31-60
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 4, 4, & ! 61-90
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 0, 0/) ! 91-106
!
! Slater exponents for valence states
! Aleksander Herman
! Empirically adjusted and consistent set of EHT valence orbital parameters for all elements of the periodic table
! Modelling and Simulation in Materials Science and Engineering, 12, 21-32 (2004)
!
! Hydrogen uses 1.2000, not the original 1.0000
!
REAL(KIND=dp), DIMENSION(4, 106) :: sexp = RESHAPE((/ &
1.2000, 0.0000, 0.0000, 0.0000, 1.6469, 0.0000, 0.0000, 0.0000, 0.6534, 0.5305, 0.0000, 0.0000, & ! 1 2 3
1.0365, 0.8994, 0.0000, 0.0000, 1.3990, 1.2685, 0.0000, 0.0000, 1.7210, 1.6105, 0.0000, 0.0000, & ! 4 5 6
2.0348, 1.9398, 0.0000, 0.0000, 2.2399, 2.0477, 0.0000, 0.0000, 2.5644, 2.4022, 0.0000, 0.0000, & ! 7 8 9
2.8812, 2.7421, 0.0000, 0.0000, 0.8675, 0.6148, 0.0000, 0.0000, 1.1935, 0.8809, 0.0000, 0.0000, & ! 10 11 12
1.5143, 1.1660, 0.0000, 0.0000, 1.7580, 1.4337, 0.0000, 0.0000, 1.9860, 1.6755, 0.0000, 0.0000, & ! 13 14 15
2.1362, 1.7721, 0.0000, 0.0000, 2.3617, 2.0176, 0.0000, 0.0000, 2.5796, 2.2501, 0.0000, 0.0000, & ! 16 17 18
0.9362, 0.6914, 0.0000, 0.0000, 1.2112, 0.9329, 0.0000, 0.0000, 1.2870, 0.9828, 2.4341, 0.0000, & ! 19 20 21
1.3416, 1.0104, 2.6439, 0.0000, 1.3570, 0.9947, 2.7809, 0.0000, 1.3804, 0.9784, 2.9775, 0.0000, & ! 22 23 24
1.4761, 1.0641, 3.2208, 0.0000, 1.5465, 1.1114, 3.4537, 0.0000, 1.5650, 1.1001, 3.6023, 0.0000, & ! 25 26 27
1.5532, 1.0594, 3.7017, 0.0000, 1.5791, 1.0527, 3.8962, 0.0000, 1.7778, 1.2448, 0.0000, 0.0000, & ! 28 29 30
2.0675, 1.5073, 0.0000, 0.0000, 2.2702, 1.7680, 0.0000, 0.0000, 2.4546, 1.9819, 0.0000, 0.0000, & ! 31 32 33
2.5680, 2.0548, 0.0000, 0.0000, 2.7523, 2.2652, 0.0000, 0.0000, 2.9299, 2.4617, 0.0000, 0.0000, & ! 34 35 36
1.0963, 0.7990, 0.0000, 0.0000, 1.3664, 1.0415, 0.0000, 0.0000, 1.4613, 1.1100, 2.1576, 0.0000, & ! 37 38 39
1.5393, 1.1647, 2.3831, 0.0000, 1.5926, 1.1738, 2.6256, 0.0000, 1.6579, 1.2186, 2.8241, 0.0000, & ! 40 41 42
1.6930, 1.2490, 2.9340, 0.0000, 1.7347, 1.2514, 3.1524, 0.0000, 1.7671, 1.2623, 3.3113, 0.0000, & ! 43 44 45
1.6261, 1.1221, 3.0858, 0.0000, 1.8184, 1.2719, 3.6171, 0.0000, 1.9900, 1.4596, 0.0000, 0.0000, & ! 46 47 48
2.4649, 1.6848, 0.0000, 0.0000, 2.4041, 1.9128, 0.0000, 0.0000, 2.5492, 2.0781, 0.0000, 0.0000, & ! 49 50 51
2.6576, 2.1718, 0.0000, 0.0000, 2.8080, 2.3390, 0.0000, 0.0000, 2.9595, 2.5074, 0.0000, 0.0000, & ! 52 53 54
1.1993, 0.8918, 0.0000, 0.0000, 1.4519, 1.1397, 0.0000, 0.0000, 1.5331, 1.1979, 2.2743, 4.4161, & ! 55 56 57
1.5379, 1.1930, 2.2912, 4.9478, 1.5162, 1.1834, 2.0558, 4.8982, 1.5322, 1.1923, 2.0718, 5.0744, & ! 58 59 60
1.5486, 1.2018, 2.0863, 5.2466, 1.5653, 1.2118, 2.0999, 5.4145, 1.5762, 1.2152, 2.0980, 5.5679, & ! 61 62 63
1.6703, 1.2874, 2.4862, 5.9888, 1.6186, 1.2460, 2.1383, 5.9040, 1.6358, 1.2570, 2.1472, 6.0598, & ! 64 65 66
1.6536, 1.2687, 2.1566, 6.2155, 1.6723, 1.2813, 2.1668, 6.3703, 1.6898, 1.2928, 2.1731, 6.5208, & ! 67 68 69
1.7063, 1.3030, 2.1754, 6.6686, 1.6647, 1.2167, 2.3795, 0.0000, 1.8411, 1.3822, 2.7702, 0.0000, & ! 70 71 72
1.9554, 1.4857, 3.0193, 0.0000, 2.0190, 1.5296, 3.1936, 0.0000, 2.0447, 1.5276, 3.3237, 0.0000, & ! 73 74 75
2.1361, 1.6102, 3.5241, 0.0000, 2.2167, 1.6814, 3.7077, 0.0000, 2.2646, 1.6759, 3.8996, 0.0000, & ! 76 77 78
2.3185, 1.7126, 4.0525, 0.0000, 2.4306, 1.8672, 0.0000, 0.0000, 2.5779, 1.9899, 0.0000, 0.0000, & ! 79 80 81
2.7241, 2.1837, 0.0000, 0.0000, 2.7869, 2.2146, 0.0000, 0.0000, 2.9312, 2.3830, 0.0000, 0.0000, & ! 82 83 84
3.1160, 2.6200, 0.0000, 0.0000, 3.2053, 2.6866, 0.0000, 0.0000, 1.4160, 1.0598, 0.0000, 0.0000, & ! 85 86 87
1.6336, 1.3011, 0.0000, 0.0000, 1.6540, 1.2890, 2.3740, 3.7960, 1.8381, 1.4726, 2.6584, 4.3613, & ! 88 89 90
1.7770, 1.4120, 2.5710, 4.5540, 1.8246, 1.4588, 2.6496, 4.7702, 1.8451, 1.4739, 2.6940, 4.9412, & ! 91 92 93
1.7983, 1.4366, 2.5123, 4.9882, 1.8011, 1.4317, 2.5170, 5.1301, 1.8408, 1.4418, 2.7349, 5.3476, & ! 94 95 96
1.8464, 1.4697, 2.5922, 5.4596, 1.8647, 1.4838, 2.6205, 5.6140, 1.8890, 1.5050, 2.6590, 5.7740, & ! 97 98 99
1.9070, 1.5190, 2.6850, 5.9220, 1.9240, 1.5320, 2.7090, 6.0690, 1.9400, 1.5440, 2.7300, 6.2130, & ! 100 101 102
2.1300, 1.7200, 2.9900, 0.0000, 1.9200, 1.4500, 2.9700, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, & ! 103 104 105
0.0000, 0.0000, 0.0000, 0.0000/), & ! 106
(/4, 106/))
PUBLIC :: qs_gcp_env_set, qs_gcp_init
! **************************************************************************************************
CONTAINS
! **************************************************************************************************
!> \brief ...
!> \param gcp_env ...
!> \param xc_section ...
! **************************************************************************************************
SUBROUTINE qs_gcp_env_set(gcp_env, xc_section)
TYPE(qs_gcp_type), POINTER :: gcp_env
TYPE(section_vals_type), POINTER :: xc_section
CHARACTER(len=*), PARAMETER :: routineN = 'qs_gcp_env_set', routineP = moduleN//':'//routineN
CHARACTER(LEN=default_string_length), &
DIMENSION(:), POINTER :: tmpstringlist
INTEGER :: i_rep, n_rep
LOGICAL :: explicit
REAL(dp), POINTER :: params(:)
TYPE(section_vals_type), POINTER :: gcp_section
CPASSERT(ASSOCIATED(gcp_env))
gcp_section => section_vals_get_subs_vals(xc_section, "GCP_POTENTIAL")
CALL section_vals_get(gcp_section, explicit=explicit)
IF (explicit) THEN
CALL section_vals_val_get(gcp_section, "VERBOSE", l_val=gcp_env%verbose)
gcp_env%do_gcp = .TRUE.
CALL section_vals_val_get(gcp_section, "PARAMETER_FILE_NAME", &
c_val=gcp_env%parameter_file_name)
CALL section_vals_val_get(gcp_section, "GLOBAL_PARAMETERS", r_vals=params)
gcp_env%sigma = params(1)
gcp_env%alpha = params(2)
gcp_env%beta = params(3)
gcp_env%eta = params(4)
! eamiss definitions
CALL section_vals_val_get(gcp_section, "DELTA_ENERGY", n_rep_val=n_rep)
IF (n_rep > 0) THEN
ALLOCATE (gcp_env%kind_type(n_rep))
ALLOCATE (gcp_env%ea(n_rep))
DO i_rep = 1, n_rep
CALL section_vals_val_get(gcp_section, "DELTA_ENERGY", i_rep_val=i_rep, &
c_vals=tmpstringlist)
READ (tmpstringlist(1), *) gcp_env%kind_type(i_rep)
READ (tmpstringlist(2), *) gcp_env%ea(i_rep)
END DO
END IF
ELSE
gcp_env%do_gcp = .FALSE.
END IF
END SUBROUTINE qs_gcp_env_set
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param gcp_env ...
! **************************************************************************************************
SUBROUTINE qs_gcp_init(qs_env, gcp_env)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(qs_gcp_type), POINTER :: gcp_env
CHARACTER(len=*), PARAMETER :: routineN = 'qs_gcp_init', routineP = moduleN//':'//routineN
REAL(KIND=dp), PARAMETER :: epsc = 1.e-6_dp
CHARACTER(LEN=10) :: aname
CHARACTER(LEN=2) :: element_symbol
INTEGER :: i, ikind, nbas, nel, nkind, nsto, za
LOGICAL :: at_end
REAL(KIND=dp) :: ea
REAL(KIND=dp), DIMENSION(10) :: al, cl
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(cp_parser_type), POINTER :: parser
TYPE(gto_basis_set_type), POINTER :: orb_basis
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_kind_type), POINTER :: qs_kind
IF (gcp_env%do_gcp) THEN
CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, nkind=nkind)
ALLOCATE (gcp_env%gcp_kind(nkind))
DO ikind = 1, nkind
qs_kind => qs_kind_set(ikind)
gcp_env%gcp_kind(ikind)%rcsto = 0.0_dp
CALL get_qs_kind(qs_kind, element_symbol=element_symbol)
CALL get_ptable_info(element_symbol, number=za)
gcp_env%gcp_kind(ikind)%za = za
gcp_env%gcp_kind(ikind)%asto = gcp_env%eta*SUM(sexp(1:4, za))/REAL(nll(za), KIND=dp)
gcp_env%gcp_kind(ikind)%nq = nshell(za)
gcp_env%gcp_kind(ikind)%rcsto = ((nshell(za)-1)*2.5_dp-LOG(epsc))/gcp_env%gcp_kind(ikind)%asto
! basis
NULLIFY (orb_basis)
CALL get_qs_kind(qs_kind, basis_set=orb_basis, basis_type="ORB")
CALL get_gto_basis_set(gto_basis_set=orb_basis, nsgf=nbas)
nel = SUM(qs_kind%elec_conf)
gcp_env%gcp_kind(ikind)%nbvirt = REAL(nbas, KIND=dp)-0.5_dp*REAL(nel, KIND=dp)
! STO-nG
nsto = SIZE(gcp_env%gcp_kind(ikind)%al)
CALL get_sto_ng(gcp_env%gcp_kind(ikind)%asto, nsto, nshell(za), 0, al, cl)
DO i = 1, nsto
gcp_env%gcp_kind(ikind)%al(i) = al(i)
gcp_env%gcp_kind(ikind)%cl(i) = cl(i)*(2._dp*al(i)/pi)**0.75_dp
END DO
END DO
! eamiss from data file
IF (gcp_env%parameter_file_name /= "---") THEN
CALL get_qs_env(qs_env, para_env=para_env)
DO ikind = 1, nkind
qs_kind => qs_kind_set(ikind)
CALL get_qs_kind(qs_kind, element_symbol=element_symbol)
CALL get_ptable_info(element_symbol, number=za)
!
NULLIFY (parser)
CALL parser_create(parser, gcp_env%parameter_file_name, para_env=para_env)
ea = 0.0_dp
DO
at_end = .FALSE.
CALL parser_get_next_line(parser, 1, at_end)
IF (at_end) EXIT
CALL parser_get_object(parser, aname)
IF (TRIM(aname) == element_symbol) THEN
CALL parser_get_object(parser, ea)
EXIT
END IF
END DO
CALL parser_release(parser)
gcp_env%gcp_kind(ikind)%eamiss = ea
END DO
END IF
!
! eamiss from input
IF (ASSOCIATED(gcp_env%kind_type)) THEN
DO i = 1, SIZE(gcp_env%kind_type)
IF (TRIM(gcp_env%kind_type(i)) == "XX") CYCLE
element_symbol = TRIM(gcp_env%kind_type(i))
CALL get_ptable_info(element_symbol, number=za)
ea = gcp_env%ea(i)
DO ikind = 1, nkind
IF (za == gcp_env%gcp_kind(ikind)%za) THEN
gcp_env%gcp_kind(ikind)%eamiss = ea
END IF
END DO
END DO
END IF
END IF
END SUBROUTINE qs_gcp_init
! **************************************************************************************************
END MODULE qs_gcp_utils

2
src/qs_ks_methods.F
View file

@ -842,7 +842,7 @@ CONTAINS
! Sum all energy terms to obtain the total energy
energy%total = energy%core_overlap+energy%core_self+energy%core+energy%hartree+ &
energy%hartree_1c+energy%exc+energy%exc1+energy%ex+ &
energy%dispersion+energy%qmmm_el+energy%mulliken+ &
energy%dispersion+energy%gcp+energy%qmmm_el+energy%mulliken+ &
SUM(energy%ddapc_restraint)+energy%s2_restraint+ &
energy%dft_plus_u+energy%kTS+ &
energy%efield+energy%efield_core+energy%ee+ &

3
src/qs_ks_utils.F
View file

@ -1061,6 +1061,9 @@ CONTAINS
IF (energy%dispersion /= 0.0_dp) &
WRITE (UNIT=output_unit, FMT="(T3,A,T61,F20.10)") &
"Dispersion energy: ", energy%dispersion
IF (energy%gcp /= 0.0_dp) &
WRITE (UNIT=output_unit, FMT="(T3,A,T61,F20.10)") &
"gCP energy: ", energy%gcp
IF (dft_control%qs_control%gapw) THEN
WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") &
"GAPW| Exc from hard and soft atomic rho1: ", energy%exc1, &

26
src/qs_neighbor_lists.F
View file

@ -73,6 +73,7 @@ MODULE qs_neighbor_lists
USE qs_dispersion_types, ONLY: qs_dispersion_type
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_gcp_types, ONLY: qs_gcp_type
USE qs_kind_types, ONLY: get_qs_kind,&
get_qs_kind_set,&
qs_kind_type
@ -323,13 +324,14 @@ CONTAINS
TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d
TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
TYPE(neighbor_list_set_p_type), DIMENSION(:), POINTER :: saa_list, sab_all, sab_almo, &
sab_aux_fit, sab_aux_fit_asymm, sab_aux_fit_vs_orb, sab_cn, sab_core, sab_kp, sab_lrc, &
sab_orb, sab_scp, sab_se, sab_tbe, sab_vdw, sab_xb, sac_ae, sac_lri, sac_ppl, sap_oce, &
sap_ppnl, soa_list, soo_list
sab_aux_fit, sab_aux_fit_asymm, sab_aux_fit_vs_orb, sab_cn, sab_core, sab_gcp, sab_kp, &
sab_lrc, sab_orb, sab_scp, sab_se, sab_tbe, sab_vdw, sab_xb, sac_ae, sac_lri, sac_ppl, &
sap_oce, sap_ppnl, soa_list, soo_list
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(paw_proj_set_type), POINTER :: paw_proj
TYPE(qs_dftb_atom_type), POINTER :: dftb_atom
TYPE(qs_dispersion_type), POINTER :: dispersion_env
TYPE(qs_gcp_type), POINTER :: gcp_env
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(section_vals_type), POINTER :: hfx_sections, neighbor_list_section
@ -827,6 +829,24 @@ CONTAINS
END IF
END IF
! Build the neighbor lists for the gCP pair potential
NULLIFY (gcp_env)
CALL get_qs_env(qs_env=qs_env, gcp_env=gcp_env)
IF (ASSOCIATED(gcp_env)) THEN
IF (gcp_env%do_gcp) THEN
sab_gcp => gcp_env%sab_gcp
DO ikind = 1, nkind
c_radius(ikind) = gcp_env%gcp_kind(ikind)%rcsto
END DO
CALL pair_radius_setup(orb_present, orb_present, c_radius, c_radius, pair_radius)
CALL build_neighbor_lists(sab_gcp, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, operator_type="PP", nlname="sab_gcp")
gcp_env%sab_gcp => sab_gcp
ELSE
NULLIFY (gcp_env%sab_gcp)
END IF
END IF
IF (lrigpw .OR. lri_optbas) THEN
! set neighborlists in lri_env environment
CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius)

3
src/qs_scf_output.F
View file

@ -425,6 +425,9 @@ CONTAINS
IF (energy%dispersion /= 0.0_dp) &
WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14)") &
"Dispersion energy: ", energy%dispersion
IF (energy%gcp /= 0.0_dp) &
WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14)") &
"gCP energy: ", energy%gcp
IF (gapw) THEN
WRITE (UNIT=output_unit, FMT="(/,(T3,A,T56,F25.14))") &
"GAPW| Exc from hard and soft atomic rho1: ", energy%exc1, &

5
src/subsys/atprop_types.F
View file

@ -34,6 +34,7 @@ MODULE atprop_types
REAL(KIND=dp), DIMENSION(:), POINTER :: atexc
REAL(KIND=dp), DIMENSION(:), POINTER :: ateself
REAL(KIND=dp), DIMENSION(:), POINTER :: atevdw
REAL(KIND=dp), DIMENSION(:), POINTER :: ategcp
REAL(KIND=dp), DIMENSION(:), POINTER :: atecc
REAL(KIND=dp), DIMENSION(:), POINTER :: ate1c
REAL(KIND=dp), DIMENSION(:), POINTER :: atecoul
@ -54,7 +55,7 @@ CONTAINS
CALL atprop_release(atprop_env)
ALLOCATE (atprop_env)
NULLIFY (atprop_env%atener, atprop_env%atstress)
NULLIFY (atprop_env%ateb, atprop_env%atevdw, atprop_env%atecc, atprop_env%atecoul)
NULLIFY (atprop_env%ateb, atprop_env%atevdw, atprop_env%ategcp, atprop_env%atecc, atprop_env%atecoul)
NULLIFY (atprop_env%ateself, atprop_env%atexc, atprop_env%ate1c)
atprop_env%energy = .FALSE.
atprop_env%stress = .FALSE.
@ -78,6 +79,7 @@ CONTAINS
CALL atprop_array_init(atprop_env%atener, natom)
CALL atprop_array_release(atprop_env%ateb)
CALL atprop_array_release(atprop_env%atevdw)
CALL atprop_array_release(atprop_env%ategcp)
CALL atprop_array_release(atprop_env%atecc)
CALL atprop_array_release(atprop_env%atecoul)
CALL atprop_array_release(atprop_env%ateself)
@ -168,6 +170,7 @@ CONTAINS
CALL atprop_array_release(atprop_env%ateself)
CALL atprop_array_release(atprop_env%atexc)
CALL atprop_array_release(atprop_env%atevdw)
CALL atprop_array_release(atprop_env%ategcp)
CALL atprop_array_release(atprop_env%atecc)
CALL atprop_array_release(atprop_env%ate1c)
CALL atprop_array_release(atprop_env%atecoul)

99
tests/QS/regtest-hf-3c/ch2o.inp Normal file
View file

@ -0,0 +1,99 @@
&GLOBAL
PROJECT ch2o
RUN_TYPE ENERGY
PRINT_LEVEL LOW
&END GLOBAL
&FORCE_EVAL
METHOD QS
&DFT
BASIS_SET_FILE_NAME BASIS_MINIX
POTENTIAL_FILE_NAME POTENTIAL
&MGRID
CUTOFF 200
REL_CUTOFF 30
&END MGRID
&QS
METHOD GAPW
&END QS
&SCF
SCF_GUESS ATOMIC
MAX_SCF 1
EPS_SCF 1.0e-6
&END SCF
&XC
&XC_FUNCTIONAL NONE
&END XC_FUNCTIONAL
&HF
&SCREENING
EPS_SCHWARZ 1.0E-10
&END
&END
&vdW_POTENTIAL
DISPERSION_FUNCTIONAL PAIR_POTENTIAL
&PAIR_POTENTIAL
TYPE DFTD3(BJ)
D3BJ_SCALING 1.0000 0.4171 0.8777 2.9149
PARAMETER_FILE_NAME dftd3.dat
SHORT_RANGE_CORRECTION TRUE
SHORT_RANGE_CORRECTION_PARAMETERS 0.03 0.70 1.50 0.75
&PRINT_DFTD
&END
&END PAIR_POTENTIAL
&END vdW_POTENTIAL
&gcp_potential
GLOBAL_PARAMETERS 0.1290 1.1549 1.1763 1.1526
DELTA_ENERGY H 0.04240
DELTA_ENERGY O 0.47901
DELTA_ENERGY C 0.27995
DELTA_ENERGY N 0.35791
VERBOSE T
# elem emiss nbas elem emiss nbas elem emiss nbas
# h 0.04240 1 he 0.02832 1 li 0.17787 5
# be 0.17160 5 b 0.22424 5 c 0.27995 5
# n 0.35791 5 o 0.47901 5 f 0.63852 5
# ne 0.83235 5 na 1.11411 9 mg 1.27115 9
# al 1.44695 14 si 1.61098 14 p 1.76661 14
# s 1.98823 14 cl 2.22845 14 ar 2.48796 14
# k 0.37425 11 ca 0.46097 11 sc 0.44489 21
# ti 0.40499 21 v 0.37841 21 cr 0.37344 21
# mn 0.36125 21 fe 0.36001 21 co 0.36293 21
# ni 0.24380 21 cu 0.40530 21 zn 0.39651 21
# ga 0.35002 32 ge 0.34578 32 as 0.34953 32
# se 0.36731 32 br 0.38201 32 kr 0.39971 32
&END GCP_POTENTIAL
&END XC
&POISSON
PERIODIC NONE
PSOLVER MT
&END
&END DFT
&SUBSYS
&CELL
ABC 10.0 10.0 10.0
PERIODIC NONE
&END CELL
&COORD
UNIT bohr
O 0.097071 0.000000 0.000000
C 2.416225 0.000000 0.000000
H 3.534655 1.775600 0.000000
H 3.534655 -1.775600 0.000000
&END COORD
&KIND H
BASIS_SET minix
POTENTIAL ALL
&END
&KIND C
BASIS_SET minix
POTENTIAL ALL
&END
&KIND O
BASIS_SET minix
POTENTIAL ALL
&END
&KIND N
BASIS_SET minix
POTENTIAL ALL
&END
&END SUBSYS
&END FORCE_EVAL

115
tests/QS/regtest-hf-3c/exam.inp Normal file
View file

@ -0,0 +1,115 @@
&GLOBAL
PROJECT s22-12
RUN_TYPE ENERGY
PRINT_LEVEL LOW
&END GLOBAL
&FORCE_EVAL
METHOD QS
&DFT
BASIS_SET_FILE_NAME BASIS_MINIX
POTENTIAL_FILE_NAME POTENTIAL
&MGRID
CUTOFF 200
REL_CUTOFF 30
&END MGRID
&QS
METHOD GAPW
&END QS
&SCF
SCF_GUESS ATOMIC
MAX_SCF 1
EPS_SCF 1.0e-6
&END SCF
&XC
&XC_FUNCTIONAL NONE
&END XC_FUNCTIONAL
&HF
&SCREENING
EPS_SCHWARZ 1.0E-10
&END
&END
&vdW_POTENTIAL
DISPERSION_FUNCTIONAL PAIR_POTENTIAL
&PAIR_POTENTIAL
TYPE DFTD3(BJ)
D3BJ_SCALING 1.0000 0.4171 0.8777 2.9149
PARAMETER_FILE_NAME dftd3.dat
SHORT_RANGE_CORRECTION TRUE
SHORT_RANGE_CORRECTION_PARAMETERS 0.03 0.70 1.50 0.75
&PRINT_DFTD
&END
&END PAIR_POTENTIAL
&END vdW_POTENTIAL
&gcp_potential
GLOBAL_PARAMETERS 0.1290 1.1549 1.1763 1.1526
DELTA_ENERGY H 0.04240
DELTA_ENERGY O 0.47901
DELTA_ENERGY C 0.27995
DELTA_ENERGY N 0.35791
VERBOSE T
# elem emiss nbas elem emiss nbas elem emiss nbas
# h 0.04240 1 he 0.02832 1 li 0.17787 5
# be 0.17160 5 b 0.22424 5 c 0.27995 5
# n 0.35791 5 o 0.47901 5 f 0.63852 5
# ne 0.83235 5 na 1.11411 9 mg 1.27115 9
# al 1.44695 14 si 1.61098 14 p 1.76661 14
# s 1.98823 14 cl 2.22845 14 ar 2.48796 14
# k 0.37425 11 ca 0.46097 11 sc 0.44489 21
# ti 0.40499 21 v 0.37841 21 cr 0.37344 21
# mn 0.36125 21 fe 0.36001 21 co 0.36293 21
# ni 0.24380 21 cu 0.40530 21 zn 0.39651 21
# ga 0.35002 32 ge 0.34578 32 as 0.34953 32
# se 0.36731 32 br 0.38201 32 kr 0.39971 32
&END GCP_POTENTIAL
&END XC
&POISSON
PERIODIC NONE
PSOLVER MT
&END
&END DFT
&SUBSYS
&CELL
ABC 10.0 10.0 10.0
PERIODIC NONE
&END CELL
&COORD
UNIT bohr
c -2.41190081140813 -2.25210745025525 -1.32102571135184 c
c -2.41190076023348 -2.25210734627092 1.32102587176460 c
c 1.39709406754529 -4.26618935000003 1.32516272888583 c
c 1.39709422665700 -4.26618905552235 -1.32516255705396 c
c -0.61999300132123 3.95356682694306 2.15436183806086 c
c 1.63540872567219 2.56929858230983 2.15433467128730 c
c 1.63540857286923 2.56929845982578 -2.15433486995303 c
c -0.61999313131748 3.95356673906409 -2.15436197357171 c
h 2.95977650901829 -5.08196917972001 2.37710806473911 h
h 2.95977833118661 -5.08196571483010 -2.37710788056673 h
h -3.96620937903017 -1.42352690648999 -2.37446790169673 h
h -3.96620936093822 -1.42352686237886 2.37446806111936 h
h -1.52029181309796 4.49846977023347 3.91642852225071 h
h 2.52473073150291 2.00438346465477 3.91667000243445 h
h 2.52473044350564 2.00438323910073 -3.91667022666881 h
h -1.52029202228550 4.49846964283464 -3.91642862967092 h
n 2.75335234162936 1.88285287748725 -0.00000011871148 n
n -0.50612492953263 -3.26275552083227 2.63481105503609 n
n -0.50612518600038 -3.26275602762472 -2.63481089551030 n
n -1.73833355442148 4.63880381147098 -0.00000005082273 n
&END COORD
&KIND H
BASIS_SET minix
POTENTIAL ALL
&END
&KIND C
BASIS_SET minix
POTENTIAL ALL
&END
&KIND O
BASIS_SET minix
POTENTIAL ALL
&END
&KIND N
BASIS_SET minix
POTENTIAL ALL
&END
&END SUBSYS
&END FORCE_EVAL

90
tests/QS/regtest-hf-3c/h2o.inp Normal file
View file

@ -0,0 +1,90 @@
&GLOBAL
PROJECT h2o
RUN_TYPE ENERGY
PRINT_LEVEL LOW
&END GLOBAL
&FORCE_EVAL
METHOD QS
&DFT
BASIS_SET_FILE_NAME BASIS_MINIX
POTENTIAL_FILE_NAME POTENTIAL
&MGRID
CUTOFF 200
REL_CUTOFF 30
&END MGRID
&QS
METHOD GAPW
&END QS
&SCF
SCF_GUESS ATOMIC
MAX_SCF 1
EPS_SCF 1.0e-6
&END SCF
&XC
&XC_FUNCTIONAL NONE
&END XC_FUNCTIONAL
&HF
&SCREENING
EPS_SCHWARZ 1.0E-10
&END
&MEMORY
MAX_MEMORY 10
&END
&END
&vdW_POTENTIAL
DISPERSION_FUNCTIONAL PAIR_POTENTIAL
&PAIR_POTENTIAL
TYPE DFTD3(BJ)
D3BJ_SCALING 1.0000 0.4171 0.8777 2.9149
PARAMETER_FILE_NAME dftd3.dat
SHORT_RANGE_CORRECTION TRUE
SHORT_RANGE_CORRECTION_PARAMETERS 0.03 0.70 1.50 0.75
&PRINT_DFTD
&END
&END PAIR_POTENTIAL
&END vdW_POTENTIAL
&gcp_potential
GLOBAL_PARAMETERS 0.1290 1.1549 1.1763 1.1526
DELTA_ENERGY H 0.04240
DELTA_ENERGY O 0.47901
VERBOSE T
# elem emiss nbas elem emiss nbas elem emiss nbas
# h 0.04240 1 he 0.02832 1 li 0.17787 5
# be 0.17160 5 b 0.22424 5 c 0.27995 5
# n 0.35791 5 o 0.47901 5 f 0.63852 5
# ne 0.83235 5 na 1.11411 9 mg 1.27115 9
# al 1.44695 14 si 1.61098 14 p 1.76661 14
# s 1.98823 14 cl 2.22845 14 ar 2.48796 14
# k 0.37425 11 ca 0.46097 11 sc 0.44489 21
# ti 0.40499 21 v 0.37841 21 cr 0.37344 21
# mn 0.36125 21 fe 0.36001 21 co 0.36293 21
# ni 0.24380 21 cu 0.40530 21 zn 0.39651 21
# ga 0.35002 32 ge 0.34578 32 as 0.34953 32
# se 0.36731 32 br 0.38201 32 kr 0.39971 32
&END GCP_POTENTIAL
&END XC
&POISSON
PERIODIC NONE
PSOLVER MT
&END
&END DFT
&SUBSYS
&CELL
ABC 4.0 4.0 4.0
PERIODIC NONE
&END CELL
&COORD
O 0.000000 0.000000 -0.065587
H 0.000000 -0.757136 0.520545
H 0.000000 0.757136 0.520545
&END COORD
&KIND H
BASIS_SET minix
POTENTIAL ALL
&END
&KIND O
BASIS_SET minix
POTENTIAL ALL
&END
&END SUBSYS
&END FORCE_EVAL

99
tests/QS/regtest-hf-3c/n3.inp Normal file
View file

@ -0,0 +1,99 @@
&GLOBAL
PROJECT n3
RUN_TYPE ENERGY
PRINT_LEVEL LOW
&END GLOBAL
&FORCE_EVAL
METHOD QS
&DFT
LSD
BASIS_SET_FILE_NAME BASIS_MINIX
POTENTIAL_FILE_NAME POTENTIAL
&MGRID
CUTOFF 200
REL_CUTOFF 30
&END MGRID
&QS
METHOD GAPW
&END QS
&SCF
SCF_GUESS ATOMIC
MAX_SCF 1
EPS_SCF 1.0e-6
&END SCF
&XC
&XC_FUNCTIONAL NONE
&END XC_FUNCTIONAL
&HF
&SCREENING
EPS_SCHWARZ 1.0E-10
&END
&END
&vdW_POTENTIAL
DISPERSION_FUNCTIONAL PAIR_POTENTIAL
&PAIR_POTENTIAL
TYPE DFTD3(BJ)
D3BJ_SCALING 1.0000 0.4171 0.8777 2.9149
PARAMETER_FILE_NAME dftd3.dat
SHORT_RANGE_CORRECTION TRUE
SHORT_RANGE_CORRECTION_PARAMETERS 0.03 0.70 1.50 0.75
&PRINT_DFTD
&END
&END PAIR_POTENTIAL
&END vdW_POTENTIAL
&gcp_potential
GLOBAL_PARAMETERS 0.1290 1.1549 1.1763 1.1526
DELTA_ENERGY H 0.04240
DELTA_ENERGY O 0.47901
DELTA_ENERGY C 0.27995
DELTA_ENERGY N 0.35791
VERBOSE T
# elem emiss nbas elem emiss nbas elem emiss nbas
# h 0.04240 1 he 0.02832 1 li 0.17787 5
# be 0.17160 5 b 0.22424 5 c 0.27995 5
# n 0.35791 5 o 0.47901 5 f 0.63852 5
# ne 0.83235 5 na 1.11411 9 mg 1.27115 9
# al 1.44695 14 si 1.61098 14 p 1.76661 14
# s 1.98823 14 cl 2.22845 14 ar 2.48796 14
# k 0.37425 11 ca 0.46097 11 sc 0.44489 21
# ti 0.40499 21 v 0.37841 21 cr 0.37344 21
# mn 0.36125 21 fe 0.36001 21 co 0.36293 21
# ni 0.24380 21 cu 0.40530 21 zn 0.39651 21
# ga 0.35002 32 ge 0.34578 32 as 0.34953 32
# se 0.36731 32 br 0.38201 32 kr 0.39971 32
&END GCP_POTENTIAL
&END XC
&POISSON
PERIODIC NONE
PSOLVER MT
&END
&END DFT
&SUBSYS
&CELL
ABC 10.0 10.0 10.0
PERIODIC NONE
&END CELL
&COORD
UNIT bohr
O 0.000000 0.000000 0.000000
O 4.000000 0.000000 0.000000
O 8.000000 0.000000 0.000000
&END COORD
&KIND H
BASIS_SET minix
POTENTIAL ALL
&END
&KIND C
BASIS_SET minix
POTENTIAL ALL
&END
&KIND O
BASIS_SET minix
POTENTIAL ALL
&END
&KIND N
BASIS_SET minix
POTENTIAL ALL
&END
&END SUBSYS
&END FORCE_EVAL

98
tests/QS/regtest-hf-3c/o2.inp Normal file
View file

@ -0,0 +1,98 @@
&GLOBAL
PROJECT o2
RUN_TYPE ENERGY
PRINT_LEVEL LOW
&END GLOBAL
&FORCE_EVAL
METHOD QS
&DFT
LSD
BASIS_SET_FILE_NAME BASIS_MINIX
POTENTIAL_FILE_NAME POTENTIAL
&MGRID
CUTOFF 200
REL_CUTOFF 30
&END MGRID
&QS
METHOD GAPW
&END QS
&SCF
SCF_GUESS ATOMIC
MAX_SCF 1
EPS_SCF 1.0e-6
&END SCF
&XC
&XC_FUNCTIONAL NONE
&END XC_FUNCTIONAL
&HF
&SCREENING
EPS_SCHWARZ 1.0E-10
&END
&END
&vdW_POTENTIAL
DISPERSION_FUNCTIONAL PAIR_POTENTIAL
&PAIR_POTENTIAL
TYPE DFTD3(BJ)
D3BJ_SCALING 1.0000 0.4171 0.8777 2.9149
PARAMETER_FILE_NAME dftd3.dat
SHORT_RANGE_CORRECTION TRUE
SHORT_RANGE_CORRECTION_PARAMETERS 0.03 0.70 1.50 0.75
&PRINT_DFTD
&END
&END PAIR_POTENTIAL
&END vdW_POTENTIAL
&gcp_potential
GLOBAL_PARAMETERS 0.1290 1.1549 1.1763 1.1526
DELTA_ENERGY H 0.04240
DELTA_ENERGY O 0.47901
DELTA_ENERGY C 0.27995
DELTA_ENERGY N 0.35791
VERBOSE T
# elem emiss nbas elem emiss nbas elem emiss nbas
# h 0.04240 1 he 0.02832 1 li 0.17787 5
# be 0.17160 5 b 0.22424 5 c 0.27995 5
# n 0.35791 5 o 0.47901 5 f 0.63852 5
# ne 0.83235 5 na 1.11411 9 mg 1.27115 9
# al 1.44695 14 si 1.61098 14 p 1.76661 14
# s 1.98823 14 cl 2.22845 14 ar 2.48796 14
# k 0.37425 11 ca 0.46097 11 sc 0.44489 21
# ti 0.40499 21 v 0.37841 21 cr 0.37344 21
# mn 0.36125 21 fe 0.36001 21 co 0.36293 21
# ni 0.24380 21 cu 0.40530 21 zn 0.39651 21
# ga 0.35002 32 ge 0.34578 32 as 0.34953 32
# se 0.36731 32 br 0.38201 32 kr 0.39971 32
&END GCP_POTENTIAL
&END XC
&POISSON
PERIODIC NONE
PSOLVER MT
&END
&END DFT
&SUBSYS
&CELL
ABC 10.0 10.0 10.0
PERIODIC NONE
&END CELL
&COORD
UNIT bohr
O 0.000000 0.000000 0.000000
O 4.000000 0.000000 0.000000
&END COORD
&KIND H
BASIS_SET minix
POTENTIAL ALL
&END
&KIND C
BASIS_SET minix
POTENTIAL ALL
&END
&KIND O
BASIS_SET minix
POTENTIAL ALL
&END
&KIND N
BASIS_SET minix
POTENTIAL ALL
&END
&END SUBSYS
&END FORCE_EVAL

197
tests/QS/regtest-hf-3c/taxol.inp Normal file
View file

@ -0,0 +1,197 @@
&GLOBAL
PROJECT taxol
RUN_TYPE ENERGY
PRINT_LEVEL LOW
&END GLOBAL
&FORCE_EVAL
METHOD QS
&DFT
BASIS_SET_FILE_NAME BASIS_MINIX
POTENTIAL_FILE_NAME POTENTIAL
&MGRID
CUTOFF 200
REL_CUTOFF 30
&END MGRID
&QS
METHOD GAPW
EPS_DEFAULT 1.e-16
&END QS
&SCF
SCF_GUESS ATOMIC
MAX_SCF 1
EPS_SCF 1.0e-6
&END SCF
&XC
&XC_FUNCTIONAL NONE
&END XC_FUNCTIONAL
&HF
&SCREENING
EPS_SCHWARZ 1.0E-10
&END
&END
&vdW_POTENTIAL
DISPERSION_FUNCTIONAL PAIR_POTENTIAL
&PAIR_POTENTIAL
TYPE DFTD3(BJ)
D3BJ_SCALING 1.0000 0.4171 0.8777 2.9149
PARAMETER_FILE_NAME dftd3.dat
SHORT_RANGE_CORRECTION TRUE
SHORT_RANGE_CORRECTION_PARAMETERS 0.03 0.70 1.50 0.75
&PRINT_DFTD
&END
&END PAIR_POTENTIAL
&END vdW_POTENTIAL
&gcp_potential
GLOBAL_PARAMETERS 0.1290 1.1549 1.1763 1.1526
DELTA_ENERGY H 0.04240
DELTA_ENERGY O 0.47901
DELTA_ENERGY C 0.27995
DELTA_ENERGY N 0.35791
VERBOSE T
# elem emiss nbas elem emiss nbas elem emiss nbas
# h 0.04240 1 he 0.02832 1 li 0.17787 5
# be 0.17160 5 b 0.22424 5 c 0.27995 5
# n 0.35791 5 o 0.47901 5 f 0.63852 5
# ne 0.83235 5 na 1.11411 9 mg 1.27115 9
# al 1.44695 14 si 1.61098 14 p 1.76661 14
# s 1.98823 14 cl 2.22845 14 ar 2.48796 14
# k 0.37425 11 ca 0.46097 11 sc 0.44489 21
# ti 0.40499 21 v 0.37841 21 cr 0.37344 21
# mn 0.36125 21 fe 0.36001 21 co 0.36293 21
# ni 0.24380 21 cu 0.40530 21 zn 0.39651 21
# ga 0.35002 32 ge 0.34578 32 as 0.34953 32
# se 0.36731 32 br 0.38201 32 kr 0.39971 32
&END GCP_POTENTIAL
&END XC
&POISSON
PERIODIC NONE
PSOLVER MT
&END
&END DFT
&SUBSYS
&CELL
ABC 20.0 20.0 20.0
PERIODIC NONE
&END CELL
&COORD
UNIT bohr
C -2.75422175362833 4.88205595547640 -3.35668019764011
C -4.89636437694867 3.48111474498741 -4.66577349906772
C -3.14587898089591 7.00916526986790 -1.47767402736533
C -0.39112727923381 3.98833834149900 -3.64501356107012
C -4.65506670157627 0.68960657647565 -3.78685362412784
C -2.31332074784111 5.93805020424707 1.12814603019323
C 0.09089322841768 1.46235809123155 -4.91230651085006
C 1.96693005487524 5.26082579230782 -2.70107385817416
C -7.61685041311990 4.41676098862359 -4.27367586421650
O -5.62383954061628 8.03714545877456 -1.33049858660306
C -4.52045648245173 3.65946126011641 -7.55996679527227
H -1.84177110780185 8.57576662860584 -1.89303703430124
C -5.13243580268152 0.40403946264478 -0.86393188517975
C -3.71347095525440 3.79294546628795 2.51849054375829
O -6.61168374546654 -0.65386883548323 -5.01235773979409
C -2.10762359900994 -0.42898639683603 -4.65782103220549
O -0.35238413511144 6.74945463069459 1.94980804066986
O 2.22524359337752 0.24754838843429 -3.69760951132643
C -5.97657298910859 10.44092860818972 -2.08867755277763
H -5.84453199337999 2.39429290249141 -8.50026768917311
H -2.61949375929026 3.18064406703739 -8.16681889075526
H -4.90473711850306 5.58947877125332 -8.16765575917408
H -8.88162145572969 3.15570277168137 -5.30027698963047
H -7.83825477939877 6.30836465515381 -5.04920693291460
H -8.21625444059990 4.43661593457266 -2.31510552529948
H 3.48574782410364 4.97906414496096 -4.06580553999356
H 2.55751755667330 4.45198098416919 -0.89670761081202
H 1.69543966453997 7.28141481571322 -2.43557502863277
H 0.62286616389413 1.72785562695556 -6.90321258147153
C -3.03248696059022 1.28827683819439 1.03261433538677
C -2.66036475699012 3.90306380603502 5.31049897601019
C -6.55109930013046 4.33728093195026 2.78265895208614
O -5.68447896699966 -2.24003494226763 -0.49225534896911
H -6.77127058463450 -2.30742873301951 -4.24721563014199
C 4.38919515350772 -0.01727680729772 -4.92750625715577
H -2.42067549348539 -1.33700865832034 -6.48496725784320
H -1.53950069221244 -1.90354293785717 -3.33566589444156
H -6.89722109959338 1.40039398258999 -0.45980320279725
C -8.65977934424017 11.25189491756870 -1.60795349896372
O -4.36186712637548 11.73356468493966 -3.00602571174966
C -2.13846158901065 -0.93507487124965 2.71225450939085
C -0.17601450117489 2.48050180266109 5.63641123557416
O -4.46642024637622 3.09117055703745 7.08004603516665
H -1.32761907425908 1.71171379305483 -0.06004279427161
C -7.88079656414129 -2.86785289931735 0.62024273043224
H -6.79663808729591 5.91464781461330 4.08330652620012
H -7.51403163140735 2.69902445675473 3.57087654500199
H -7.44886189316021 4.84692014828003 1.01541211055408
H -2.32994502928741 5.91837360344958 5.70368188588543
C 6.23640343719723 -1.71691542074622 -3.48505893590230
O 4.90608977364205 0.96358302217710 -6.91255859364461
H -8.95018172414545 13.14957685773710 -2.33874539972598
H -9.01383089081425 11.21828386500458 0.41939820347943
H -9.95426564852871 9.93490300004304 -2.51423269320752
C -0.50442470804037 -0.30326765337150 5.04548880188861
C -3.95357943210590 -2.16405001923217 4.61989078502303
O -0.96833297587134 -2.66006294543500 0.96045744368580
H -4.06828240944650 1.37755551634343 7.62563495419667
H 1.26792608465731 3.28800786551726 4.41183663823650
H 0.44887388611826 2.66812100457173 7.59556351117105
O -9.37033865414920 -1.34693826681822 1.43115546398267
C -8.23241043866238 -5.64735838430951 0.77816295700356
C 8.46218736839208 -0.05304708887146 -2.51253314803634
O 5.02077123832797 -3.02350801140142 -1.53581458156187
H 6.99182242704328 -3.10269620933609 -4.83385300940165
O -2.33116640649002 -1.49977678243844 6.67761982092692
C 0.01884127140618 -4.88468883322427 1.64225736838665
H -4.16500289381031 -4.21140182964611 4.43716572795657
H -5.80481804486070 -1.25624563717105 4.78335998170243
H 1.30702395724853 -1.30312799473358 5.09896111525751
C -10.31931343489841 -6.56118929268801 2.10750154472034
C -6.53872425052417 -7.33640405183194 -0.33144081060646
C 10.75166785917711 -1.53835912450025 -1.65885800844686
N 7.46999111101415 1.63117940305486 -0.54579410827088
H 9.03497087444481 1.17943011266750 -4.08134444025290
H 3.72272508337098 -1.95075651250539 -0.78857749276112
O 0.26038541584440 -5.62138263220708 3.77448961105768
C 0.79020534468577 -6.42909118298574 -0.62945542330994
C -10.69898451734659 -9.14038461994881 2.33461325654239
C -6.92948455519594 -9.91464865360206 -0.09902701044444
H -11.61898838013555 -5.22445034257893 2.95781253200858
H -4.92664885367749 -6.62377821613846 -1.37394969966546
C 13.06741392749892 -0.28224600190127 -1.57185466919596
C 10.64491807061757 -4.07394950684923 -0.96743510107778
C 6.30353106105665 0.59634471844069 1.57090033564930
H 6.24477407196677 2.88198412021953 -1.34346073680208
H -0.45167574558374 -8.06746588469714 -0.73469509421189
H 2.71612645607674 -7.09512152613059 -0.35626910730471
H 0.66378643268868 -5.34375908969597 -2.36648836149020
C -9.00331942114607 -10.81657464414146 1.23519826465783
H -12.31610332239828 -9.85003552359730 3.37248177184792
H -5.61190905597755 -11.22746482819311 -0.95799714348069
C 15.23250461685864 -1.53138908472102 -0.78543460344008
C 12.81584681145065 -5.32131711143592 -0.18580328467309
O 4.13643431181960 1.19055530546031 2.07477806937000
C 7.85671181743914 -1.03597887652096 3.23907719586055
H 13.16708189075178 1.69320102500146 -2.11033513523406
H 8.86964807254477 -5.08894394415900 -0.99666723342523
H -9.29830568244296 -12.83574606425507 1.41784191742182
C 15.11010391360424 -4.05696073250481 -0.08711166034559
C 10.32350874633572 -0.34179122139165 3.85153710083126
C 6.79252560972846 -3.18951151680852 4.32074038536061
H 17.01806987741279 -0.52911033379028 -0.71754578030844
H 12.70545902828425 -7.29286656425878 0.36073524265411
C 11.72953198165022 -1.82482837492992 5.48955863706336
C 8.21575004808072 -4.68013477485924 5.93615318804377
H 16.79925354836453 -5.03600851477630 0.53276304207989
H 11.11963623702429 1.36747523392382 3.05057597822599
H 4.86666579164096 -3.70660320433175 3.85606162241743
C 10.68418523664758 -4.00288306661292 6.51656221602793
H 13.64856556643712 -1.28400567524983 5.96139006086769
H 7.39457824978855 -6.36991639262988 6.75199678774226
H 11.79525679743083 -5.16879168381238 7.78319070073019
&END COORD
&KIND DEFAULT
BASIS_SET minix
POTENTIAL ALL
&END
&END SUBSYS
&END FORCE_EVAL