1!
2!   Polarizable Embedding (PE) library
3!   Copyright (C) 2013-2018, 2020 The PE library developers. See the CONTRIBUTORS file
4!                                 in the top-level directory of this distribution.
5!
6!   This file is part of the PE library.
7!
8!   The PE library is free software: you can redistribute it and/or modify
9!   it under the terms of the GNU Lesser General Public License as
10!   published by the Free Software Foundation, either version 3 of the
11!   License, or (at your option) any later version.
12!
13!   The PE library is distributed in the hope that it will be useful,
14!   but WITHOUT ANY WARRANTY; without even the implied warranty of
15!   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16!   GNU Lesser General Public License for more details.
17!
18!   You should have received a copy of the GNU Lesser General Public License
19!   along with the PE library. If not, see <http://www.gnu.org/licenses/>.
20!
21!   Contact information:
22!
23!   Jogvan Magnus Haugaard Olsen
24!   E-mail: foeroyingur@gmail.com
25!
26module pelib_options
27
28    use pelib_precision
29
30    implicit none
31
32    ! options
33    logical(lp), public, save :: pelib_gspol = .false.
34    logical(lp), public, save :: pelib_iter = .true.
35    logical(lp), public, save :: pelib_diis = .false.
36    logical(lp), public, save :: pelib_redthr = .false.
37    logical(lp), public, save :: pelib_border = .false.
38    logical(lp), public, save :: pelib_ind_damp = .false.
39    logical(lp), public, save :: pelib_mul_damp = .false.
40    logical(lp), public, save :: pelib_core_damp = .false.
41    logical(lp), public, save :: pelib_amoeba_damp = .false.
42    logical(lp), public, save :: pelib_nomb = .false.
43    logical(lp), public, save :: pelib_polar = .false.
44    logical(lp), public, save :: pelib_cube = .false.
45    logical(lp), public, save :: pelib_mep = .false.
46    logical(lp), public, save :: pelib_skipqm = .false.
47    logical(lp), public, save :: pelib_twoints = .false.
48    logical(lp), public, save :: pelib_repuls = .true.
49    logical(lp), public, save :: pelib_savden = .false.
50    logical(lp), public, save :: pelib_fd = .false.
51    logical(lp), public, save :: pelib_fdes = .true.
52    logical(lp), public, save :: pelib_sol = .false.
53    logical(lp), public, save :: pelib_noneq = .true.
54    logical(lp), public, save :: pelib_restart = .false.
55    logical(lp), public, save :: pelib_verbose = .false.
56    logical(lp), public, save :: pelib_debug = .false.
57    logical(lp), public, save :: pelib_fixsol = .false.
58    logical(lp), public, save :: pelib_nomuq = .false.
59    logical(lp), public, save :: pelib_novmu = .false.
60    logical(lp), public, save :: pelib_read_surf = .false.
61    logical(lp), public, save :: pelib_zeromul = .false.
62    logical(lp), public, save :: pelib_zeropol = .false.
63    logical(lp), public, save :: pelib_isopol = .false.
64    logical(lp), public, save :: pelib_gauge = .false.
65    logical(lp), public, save :: pelib_lf = .false.
66    logical(lp), public, save :: pelib_field = .false.
67    logical(lp), public, save :: pelib_skipmul = .false.
68
69    ! runtypes
70    logical(lp) :: fock = .false.
71    logical(lp) :: mixed = .false.
72    logical(lp) :: energy = .false.
73    logical(lp) :: response = .false.
74    logical(lp) :: molgrad = .false.
75    logical(lp) :: mep = .false.
76    logical(lp) :: london = .false.
77    logical(lp) :: effdipole = .false.
78    logical(lp) :: cube = .false.
79
80    ! in- and output logicals, i.e. denmats, fckmats and expvals
81    logical(lp) :: lden = .false.
82    logical(lp) :: lfck = .false.
83    logical(lp) :: lexp = .false.
84
85    ! matrix type
86    logical(lp) :: trimat = .true.
87
88    ! filenames
89    character(len=240) :: potfile = 'POTENTIAL.INP'
90    character(len=240) :: surfile = 'SURFACE.INP'
91    character(len=240) :: h5pdefile = 'standard.h5'
92    character(len=240) :: h5gridfile = 'grid.h5'
93
94    logical(lp), save :: synced = .false.
95    logical(lp), save :: initialized = .false.
96    integer(ip), save :: site_start, site_finish
97    integer(ip), save :: surp_start, surp_finish
98    integer(ip), save :: cube_start, cube_finish
99    integer(ip), dimension(:), save, allocatable :: siteloops
100    integer(ip), dimension(:), save, allocatable :: surploops
101    integer(ip), dimension(:), save, allocatable :: cubeloops
102    integer(ip), dimension(:), save, allocatable :: poldists
103    integer(ip), dimension(:), save, allocatable :: sitedists
104    integer(ip), dimension(:), save, allocatable :: surpdists
105    integer(ip), dimension(:), save, allocatable :: cubedists
106    integer(ip), dimension(:), save, allocatable :: displs
107
108    ! logical unit for output file (default is stdout)
109    integer(ip), save :: luout = 6
110
111    ! integer used to check memory allocation status,
112    ! i.e. allocate(array(len), stat=mem_stat)
113    integer(ip) :: mem_stat
114
115    ! dummy variables
116    integer(ip) :: dummy_int
117    real(rp) :: dummy_real
118    integer(ip), dimension(1) :: dummy_int_array
119    real(rp), dimension(1) :: dummy_real_array
120
121    ! constants, thresholds and stuff
122    real(rp), parameter :: zero = 1.0e-8
123    integer(ip), save :: scf_cycle = 0
124    integer(ip), save :: redlvl = 0
125    real(rp), save :: thriter = 1.0e-8
126    real(rp), save :: thrdiis = 1.0e-8
127    real(rp), save :: ind_damp = 2.1304
128    real(rp), save :: mul_damp = 2.1304
129    real(rp), save :: core_damp = 2.1304
130    real(rp), save :: amoeba_damp = 0.39
131    real(rp), save :: gauss_factor = 1.0
132    real(rp), dimension(3), save :: repfacs = 1.0
133    real(rp), save :: Rmin = 2.2
134    integer(ip), save :: nredist = 1
135    integer(ip), save :: redist_order = 1
136    character(len=6), save :: border_type = 'REDIST'
137    ! use Cholesky factorization of classical response matrix
138    logical(lp), save :: chol = .true.
139    ! solvent and dielectric constant (defaults to water)
140    character(len=80) :: solvent
141    real(rp), save :: eps = 0.0
142    real(rp), save :: epsinf = 0.0
143    ! the order from which multipoles are zeroed
144    integer(ip), save :: zeromul_order = 1
145
146    ! polarizable embedding potential info
147    ! ------------------------------------
148    ! total number of classical sites
149    integer(ip), save :: nsites = 0
150    ! number of polarizable sites
151    integer(ip), save :: npols = 0
152    ! number of surface points
153    integer(ip), save :: nsurp = 0
154    ! number fragment densities
155    integer(ip), save :: nfds = 0
156    ! number of nuclei in fragment density
157    integer(ip), save :: fdnucs = 0
158    ! multipole exclusion list length
159    integer(ip), save :: mul_lexlst = 0
160    ! polarizability exclusion list length
161    integer(ip), save :: pol_lexlst = 0
162    ! number of density matrices
163    integer(ip) :: ndens = 0
164    ! number of basis functions in core fragment
165    integer(ip), save :: nbas
166    ! size of packed matrices
167    integer(ip), save :: nnbas
168    ! size of full matrices
169    integer(ip), save :: n2bas
170    ! number of nuclei in core region
171    integer(ip), save :: nnucs = 0
172    ! number of LJ sites in core region
173    integer(ip), save :: qmLJsites = 0
174    ! index of electronic state
175    integer(ip), save :: state_idx = -1
176
177    ! Thole qm damping stuff
178    ! ----------------------
179    integer(ip), save :: nqmpols = 0
180    real(rp), dimension(:,:), allocatable, save :: qmpols
181    real(rp), save :: qmdamp = 2.1304
182
183    ! specifies what type of parameters are present
184    ! lmul(0): monopoles, lmul(1): dipoles etc.
185    logical(lp), dimension(0:5), save :: lmul = .false.
186    ! lpol(1): dipole-dipole polarizabilities
187    logical(lp), dimension(0:2,0:2), save :: lpol = .false.
188    ! lhypol(1): dipole-dipole-dipole polarizabilities/1st hyperpolarizability
189!    logical(lp), dimension(1), save :: lhypol
190!    ! lvdw: LJ parameters
191    logical(lp), save :: lvdw = .false.
192
193    ! charges, areas, coordinates, elements and exclusion lists
194    ! site nuclear charges
195    real(rp), dimension(:,:), allocatable, save :: Zs
196    ! fragment density nuclear charges
197    real(rp), dimension(:,:), allocatable, save :: Zfd
198    ! core fragment nuclear charges
199    real(rp), dimension(:,:), allocatable, save :: Zm
200    ! surface point areas
201    real(rp), dimension(:), allocatable, save :: Sa
202    ! Atom ID
203    integer(ip), dimension(:), save, allocatable :: idatm
204    ! Radii of the sphere from which each tessera stems from
205    real(rp), dimension(:), allocatable, save :: Radsp
206    ! site coordinates
207    real(rp), dimension(:,:), allocatable, save :: Rs
208    ! fragment density nuclear coordinates
209    real(rp), dimension(:,:), allocatable, save :: Rfd
210    ! core fragment nuclear coordinates
211    real(rp), dimension(:,:), allocatable, save :: Rm
212    ! all nuclear coordinates
213    !real(rp), dimension(:,:), allocatable, save :: Rm_full
214    real(rp), dimension(:,:), allocatable, save :: all_coords
215    ! surface point coordinates
216    real(rp), dimension(:,:), allocatable, save :: Rsp
217    ! The normal vector to each surface point
218    real(rp), dimension(:,:), allocatable, save :: nrmsp
219    ! core center-of-mass
220    real(rp), dimension(:), allocatable, save :: core_com
221    ! core when changing gauge
222    real(rp), dimension(:), allocatable, save :: gauge_input
223    ! core polarizabilities
224    real(rp), dimension(:,:), allocatable, save :: core_alphas
225    ! site elements
226    character(len=2), dimension(:), allocatable, save :: elems
227    ! multipole exclusion list
228    integer(ip), dimension(:,:), allocatable, save :: mul_exclists
229    ! polarizability exclusion list
230    integer(ip), dimension(:,:), allocatable, save :: pol_exclists
231
232    ! energy contributions
233    ! total
234    real(rp), dimension(:), allocatable, save :: Epe
235    ! electrostatic
236    real(rp), dimension(:,:), allocatable, save :: Ees
237    ! polarization
238    real(rp), dimension(:,:), allocatable, save :: Epol
239    ! continuum solvation
240    real(rp), dimension(:,:), allocatable, save :: Esol
241    ! fragment density
242    real(rp), dimension(:,:), allocatable, save :: Efd
243    ! LJ
244    real(rp), dimension(:), allocatable, save :: Elj
245
246    ! multipole moments
247    ! monopoles, dipoles, quadrupoles, octopoles, etc.
248    real(rp), dimension(:,:), allocatable, save :: M0s, M1s, M2s, M3s, M4s, M5s
249    ! polarizabilities
250    ! symmetric dipole-dipole polarizabilities
251    real(rp), dimension(:,:), allocatable, save :: P11s
252    ! .true. if P11 > 0 else .false.
253    logical(lp), dimension(:), allocatable, save :: zeroalphas
254
255    ! LJ parameters for MM region - from pot file
256    real(rp), dimension(:,:), allocatable, save :: LJs
257    ! LJ paramters for QM region - from dal file
258    real(rp), dimension(:,:), allocatable, save :: qmLJs
259
260    ! CUBE stuff
261    ! ---------
262
263    ! options for MEP
264    ! create QM cubes
265    logical(lp), save :: mep_qmcube = .true.
266    ! create multipole cubes
267    logical(lp), save :: mep_mulcube = .true.
268    ! external electric field
269    logical(lp), save :: mep_extfld = .false.
270    real(rp), dimension(3), save :: extfld = huge(0.0)
271    ! calculate electric field
272    logical(lp), save :: mep_field = .false.
273    ! lf component
274    logical(lp), save :: mep_lf = .false.
275    integer(ip), save :: lf_component = huge(0)
276    ! grid file
277    logical(lp), save :: mep_input = .false.
278
279    ! options for CUBE
280    ! calculate electric field
281    logical(lp), save :: cube_field = .false.
282
283    ! general cube information
284    ! number of grid points
285    integer(ip), save :: npoints
286    ! grid points
287    real(rp), dimension(:,:), allocatable, save :: Rg
288    ! CUBE file origin and step sizes
289    real(rp), dimension(3), save :: origin, step
290    ! grid density in x, y and z direction
291    integer(ip), save :: xgrid = 6
292    integer(ip), save :: ygrid = 6
293    integer(ip), save :: zgrid = 6
294    ! number of steps in x, y and z direction
295    integer(ip), save :: xsteps
296    integer(ip), save :: ysteps
297    integer(ip), save :: zsteps
298    ! box size relative to molecule size
299    real(rp), save :: xsize = 8.0
300    real(rp), save :: ysize = 8.0
301    real(rp), save :: zsize = 8.0
302
303    ! Internal field stuff and locfld stuff
304    ! Coordinates on which potential and field are calculated
305    real(rp), dimension(:,:), allocatable, save :: crds
306    ! Number of coordinates (length of crds) / 3
307    integer(ip), save :: ncrds = 0
308
309    ! External field stuff
310    ! Field strengths
311    real(rp), dimension(3), save :: efields = 0.0
312
313end module pelib_options
314