1!
2! PCMSolver, an API for the Polarizable Continuum Model
3! Copyright (C) 2020 Roberto Di Remigio, Luca Frediani and contributors.
4!
5! This file is part of PCMSolver.
6!
7! PCMSolver is free software: you can redistribute it and/or modify
8! it under the terms of the GNU Lesser General Public License as published by
9! the Free Software Foundation, either version 3 of the License, or
10! (at your option) any later version.
11!
12! PCMSolver is distributed in the hope that it will be useful,
13! but WITHOUT ANY WARRANTY; without even the implied warranty of
14! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15! GNU Lesser General Public License for more details.
16!
17! You should have received a copy of the GNU Lesser General Public License
18! along with PCMSolver.  If not, see <http://www.gnu.org/licenses/>.
19!
20! For information on the complete list of contributors to the
21! PCMSolver API, see: <http://pcmsolver.readthedocs.io/>
22!
23
24!
25!     simple input reader for cavity generator
26!     written by Krzysztof Mozgawa, 2010
27!
28!     RDR, 280114. Put things in makecav.F inside here directly.
29!
30subroutine generatecavity_cpp(maxts_, maxsph_, maxvert_,         &
31    xtscor_, ytscor_, ztscor_, ar_,                              &
32    xsphcor_, ysphcor_, zsphcor_, rsph_,                         &
33    nts_, ntsirr_, nesfp_, addsph_,                              &
34    xe_, ye_, ze_, rin_, masses_, avgArea_, rsolv_, ret_,        &
35    nr_gen_, gen1_, gen2_, gen3_,                                &
36    nvert_, vert_, centr_, isphe_, pedra_, len_pedra_)           &
37    bind(C, name='generatecavity_cpp')
38
39use, intrinsic :: iso_c_binding
40use pedra_precision
41use pedra_symmetry, only: get_point_group, point_group
42use pedra_cavity, only: polyhedra_driver
43use strings, only: carray_to_fstring
44
45implicit none
46
47#include "pcm_pcmdef.inc"
48#include "pcm_mxcent.inc"
49#include "pcm_pcm.inc"
50
51integer(c_int)  :: maxts_, maxsph_, maxvert_
52real(c_double)  :: xtscor_(maxts_), ytscor_(maxts_), ztscor_(maxts_)
53real(c_double)  :: xsphcor_(maxts_), ysphcor_(maxts_), zsphcor_(maxts_), rsph_(maxts_)
54real(c_double)  :: ar_(maxts_), xe_(maxts_), ye_(maxts_), ze_(maxts_), rin_(maxts_)
55real(c_double)  :: masses_(maxts_)
56real(c_double)  :: avgArea_, rsolv_, ret_
57integer(c_int)  :: nts_, ntsirr_, nesfp_, addsph_
58integer(c_int)  :: nr_gen_, gen1_, gen2_, gen3_
59integer(c_int)  :: nvert_(maxts_)
60real(c_double)  :: vert_(maxts_ * 30), centr_(maxts_ * 30)
61integer(c_int)  :: isphe_(maxts_)
62integer(c_int) :: len_pedra_
63character(kind=c_char, len=1), intent(in) :: pedra_(len_pedra_+1)
64
65integer(c_int)    :: i, j, k, offset
66integer(c_int)    :: error_code
67integer(kind=regint_k) :: pedra_unit
68logical           :: pedra_open, pedra_exist
69real(c_double), allocatable :: vert(:, :, :), centr(:, :, :)
70character(len=len_pedra_) :: pedra
71type(point_group) :: pgroup
72
73
74pedra = carray_to_fstring(pedra_)
75!lvpri = 121201_regint_k
76! The following INQUIRE statement returns whether the file named cavity.off is
77! connected in logical variable off_open, whether the file exists in logical
78! variable off_exist, and the unit number in integer(kind=regint_k) variable off_unit
79pedra_unit = 121201_regint_k
80inquire(file = pedra, opened = pedra_open, &
81        exist = pedra_exist)
82if (pedra_exist) then
83   open(unit = pedra_unit,                   &
84       file = pedra,       &
85       status = 'unknown',       &
86       form = 'formatted',       &
87       access = 'sequential')
88   close(unit = pedra_unit, status = 'delete')
89end if
90open(unit = pedra_unit,                      &
91    file = pedra,          &
92    status = 'new',              &
93    form = 'formatted',          &
94    access = 'sequential')
95rewind(pedra_unit)
96areats = avgArea_
97icesph = 1
98iprpcm = 3
99call get_point_group(pedra_unit, pgroup, nr_gen_, gen1_, gen2_, gen3_)
100rsolv = rsolv_
101! These parameters are fixed see one of the original GePol papers
102omega = 40.0d+00
103fro = 0.7d+00
104! ret is the minimum radius of added spheres
105ret = ret_
106! nesfp is the number of initial spheres.
107! nesf is the total number of spheres: initial + added
108nesfp = nesfp_
109do i = 1, nesfp
110   xe(i) = xe_(i)
111   ye(i) = ye_(i)
112   ze(i) = ze_(i)
113   rin(i) = rin_(i)
114   alpha(i) = 1.0d0
115end do
116! Workaround to silence compilation warning
117maxsph_ = mxsp
118maxvert_= mxver
119
120! Allocate space for the arrays containing the vertices and the centers
121! of the tesserae arcs
122allocate(vert(mxts, 10, 3))
123vert = 0.0d0
124allocate(centr(mxts, 10, 3))
125centr = 0.0d0
126
127nesf = nesfp
128
129call polyhedra_driver(pgroup, vert, centr, masses_, pedra_unit, error_code)
130
131! Common block dark magic, it will disappear one day...
132nts_ = nts
133ntsirr_ = ntsirr
134! Pass the number of added spheres back, to update the GePolCavity
135! object in the right way.
136! nesf: total number of spheres; nesfp: number of original spheres
137addsph_ = nesf - nesfp
138do i = 1, nts
139   xtscor_(i) = xtscor(i)
140   ytscor_(i) = ytscor(i)
141   ztscor_(i) = ztscor(i)
142   ar_(i) = as(i)
143   xsphcor_(i) = xe(isphe(i))
144   ysphcor_(i) = ye(isphe(i))
145   zsphcor_(i) = ze(isphe(i))
146   rsph_(i) = re(isphe(i))
147   nvert_(i) = nvert(i)
148! Fill the vert_ and centr_ arrays
149   do j = 1, nvert(i)
150      do k = 1, 3
151          offset = i + j * nts + k * nts * nvert(i)
152          vert_(offset) = vert(i, j, k)
153          centr_(offset) = centr(i, j, k)
154      end do
155   end do
156end do
157
158do i = 1, nesf
159  isphe_(i) = isphe(i)
160  xe_(i)  = xe(i)
161  ye_(i)  = ye(i)
162  ze_(i)  = ze(i)
163  rin_(i) = re(i)
164end do
165
166! Clean-up
167deallocate(vert)
168deallocate(centr)
169
170write(pedra_unit, *) "Error code is ", error_code
171write(pedra_unit, *) '<<< Done with PEDRA Fortran code >>>'
172
173close(pedra_unit)
174
175end subroutine generatecavity_cpp
176