1module pelib_cavity_generators
2
3    use pelib_precision
4
5    implicit none
6
7    private
8
9    public :: fixtes
10
11    ! Maximum number of tessera
12    integer(ip), save :: MXFFTS = 0
13    ! Number of centers
14    integer(ip), save :: NCENTS = 0
15    ! Number of tessera per centers
16    integer(ip), public, save :: NTSATM = 60
17    ! Number of tessera
18    integer(ip), public, save :: NFFTS = 0
19
20contains
21
22!------------------------------------------------------------------------------
23!
24! The following code has been implemented into the gamess program by Hui Li
25! and adapted for use in dalton(pelib)
26!
27! The tesselation is FIXPVA2
28!     Nandun M. Thellamurege and Hui Li THE JOURNAL OF CHEMICAL PHYSICS 137, 246101 (2012)
29!
30!------------------------------------------------------------------------------
31subroutine fixtes(all_centers, all_charges)
32
33    use pelib_constants
34    use pelib_options
35
36    real(rp), dimension(:,:), intent(in) :: all_centers
37    real(rp), dimension(:), intent(in) :: all_charges
38
39    real(rp), dimension(960) :: AST
40    real(rp), dimension(360) :: IDUM
41    real(rp), dimension(24) :: THEV
42    real(rp), dimension(24) :: FIV
43    real(rp), dimension(3,960) :: CDTST
44    real(rp), dimension(3,20) :: VT20
45    real(rp), dimension(6,60) :: JVT1
46    real(rp), dimension(122,3) :: CV
47    real(rp), dimension(3,40,960) :: DAIT
48    integer(ip), dimension(41,960) :: IDTMP
49    real(rp), dimension(:), allocatable :: XTS, YTS, ZTS
50    real(rp), dimension(:), allocatable :: AFIX, RFIX
51    real(rp), dimension(:,:,:), allocatable :: DAI
52    integer(ip), dimension(:), allocatable :: IDATOM
53    integer(ip), dimension(:,:), allocatable :: IDDAI
54    real(rp), dimension(:,:), allocatable :: TMPTS
55    integer(ip) :: II, I, J, INUC, NFFTS, IFFAT
56    integer(ip) :: KFFTS, ITS, JJJ
57    real(rp) :: TH, FI, CTH, STH, FIR
58    real(rp), parameter :: GOLD = 1.618033988749895
59    real(rp), parameter :: ONEGOLD = 1.0 / GOLD
60    real(rp), parameter :: SQRT13 = 0.577350269189626
61
62    EQUIVALENCE (IDUM(1),JVT1(1,1))
63
64    NCENTS = size(all_charges)
65    MXFFTS = NCENTS * NTSATM
66    ALLOCATE(XTS(MXFFTS))
67    ALLOCATE(YTS(MXFFTS))
68    ALLOCATE(ZTS(MXFFTS))
69    ALLOCATE(AFIX(MXFFTS))
70    ALLOCATE(RFIX(MXFFTS))
71    ALLOCATE(IDATOM(MXFFTS))
72    ALLOCATE(IDDAI(41,MXFFTS))
73    ALLOCATE(DAI(3,40,MXFFTS))
74    ALLOCATE(TMPTS(3,MXFFTS))
75
76    THEV  = (/0.6523581398,1.107148718,1.382085796,  &
77            &  1.759506858,2.034443936,2.489234514,  &
78            &                 0.3261790699,0.5535743589, &
79            &  0.8559571251,0.8559571251,1.017221968,&
80            &  1.229116717,1.229116717,1.433327788,  &
81            &  1.570796327,1.570796327,1.708264866,  &
82            &  1.912475937,1.912475937,2.124370686,  &
83            &  2.285635528,2.285635528,2.588018295,  &
84            &  2.815413584/)
85
86    FIV   = (/ 0.6283185307,0.0000000000,                 &
87            &  0.6283185307,0.0000000000,0.6283185307,&
88            &  0.0000000000,0.6283185307,0.0000000000,&
89            &  0.2520539002,1.004583161,0.6283185307, &
90            &  0.3293628477,0.9272742138,0.0000000000,&
91            &  0.3141592654,0.9424777961,0.6283185307,&
92            &  0.2989556830,0.9576813784,0.0000000000,&
93            &  0.3762646305,0.8803724309,0.6283188307,&
94            &  0.0000000000/)
95
96    FIR   = 1.256637061
97
98    IDUM = (/1, 6, 2, 32, 36, 37, 1, 2, 3, 33, 32, 38, 1, 3, 4, 34,  &
99     &   33, 39, 1, 4, 5, 35, 34, 40, 1, 5, 6, 36, 35, 41, 7, 2, 6, 51,&
100     &   42, 37, 8, 3, 2, 47, 43, 38, 9, 4, 3, 48, 44, 39, 10, 5, 4,   &
101     &   49, 45, 40, 11, 6, 5, 50, 46, 41, 8, 2, 12, 62, 47, 52, 9,    &
102     &   3, 13, 63, 48, 53, 10, 4, 14, 64, 49, 54, 11, 5, 15, 65, 50,  &
103     &   55, 7, 6, 16, 66, 51, 56, 7, 12, 2, 42, 57, 52, 8, 13, 3,     &
104     &   43, 58, 53, 9, 14, 4, 44, 59, 54, 10, 15, 5, 45, 60, 55, 11,  &
105     &   16, 6, 46, 61, 56, 8, 12, 18, 68, 62, 77, 9, 13, 19, 69, 63,  &
106     &   78, 10, 14, 20, 70, 64, 79, 11, 15, 21, 71, 65, 80, 7, 16,    &
107     &   17, 67, 66, 81, 7, 17, 12, 57, 67, 72, 8, 18, 13, 58, 68, 73, &
108     &   9, 19, 14, 59, 69, 74, 10, 20, 15, 60, 70, 75, 11, 21, 16,    &
109     &   61, 71, 76, 22, 12, 17, 87, 82, 72, 23, 13, 18, 88, 83, 73,   &
110     &   24, 14, 19, 89, 84, 74, 25, 15, 20, 90, 85, 75, 26, 16, 21,   &
111     &   91, 86, 76, 22, 18, 12, 82, 92, 77, 23, 19, 13, 83, 93, 78,   &
112     &   24, 20, 14, 84, 94, 79, 25, 21, 15, 85, 95, 80, 26, 17, 16,   &
113     &   86, 96, 81, 22, 17, 27, 102, 87, 97, 23, 18, 28, 103, 88, 98, &
114     &   24, 19, 29, 104, 89, 99, 25, 20, 30, 105, 90, 100, 26, 21,    &
115     &   31, 106, 91, 101, 22, 28, 18, 92, 107, 98, 23, 29, 19, 93,    &
116     &   108, 99, 24, 30, 20, 94, 109, 100, 25, 31, 21, 95, 110, 101,  &
117     &   26, 27, 17, 96, 111, 97, 22, 27, 28, 107, 102, 112, 23, 28,   &
118     &   29, 108, 103, 113, 24, 29, 30, 109, 104, 114, 25, 30, 31,     &
119     &   110, 105, 115, 26, 31, 27, 111, 106, 116, 122, 28, 27, 117,   &
120     &   118, 112, 122, 29, 28, 118, 119, 113, 122, 30, 29, 119, 120,  &
121     &   114, 122, 31, 30, 120, 121, 115, 122, 27, 31, 121, 117, 116 /)
122
123!
124!     ADAPTED FROM PCM CODE (IN GAMESS)
125!     HUI LI, NOV 26, 2011
126!
127!     COORDINATES OF VERTICES OF TESSERAE IN A SPHERE WITH UNIT RADIUS.
128!
129!                                    1
130!
131!                                 4     5
132!
133!                              3     6     2
134!
135      CV(  1,1) =  0.0
136      CV(  1,2) =  0.0
137      CV(  1,3) =  1.0
138      CV(122,1) =  0.0
139      CV(122,2) =  0.0
140      CV(122,3) = -1.0
141      II=1
142      DO I=1,24
143         TH=THEV(I)
144         FI=FIV(I)
145         CTH=COS(TH)
146         STH=SIN(TH)
147         DO J=1,5
148            FI=FI+FIR
149            IF(J.EQ.1) FI=FIV(I)
150            II=II+1
151            CV(II,1)=STH*COS(FI)
152            CV(II,2)=STH*SIN(FI)
153            CV(II,3)=CTH
154         END DO
155      END DO
156!
157      IF(NTSATM.EQ.4) THEN
158         VT20(1,1)  = 1.0
159         VT20(2,1)  = 1.0
160         VT20(3,1)  = 1.0
161         VT20(1,2)  =-1.0
162         VT20(2,2)  =-1.0
163         VT20(3,2)  = 1.0
164         VT20(1,3)  =-1.0
165         VT20(2,3)  = 1.0
166         VT20(3,3)  =-1.0
167         VT20(1,4)  = 1.0
168         VT20(2,4)  =-1.0
169         VT20(3,4)  =-1.0
170         CALL DSCAL(12,SQRT13,VT20,1)
171      END IF
172!
173      IF(NTSATM.EQ.6) THEN
174         VT20(1,1)  = 1.0
175         VT20(2,1)  = 0.0
176         VT20(3,1)  = 0.0
177         VT20(1,2)  =-1.0
178         VT20(2,2)  = 0.0
179         VT20(3,2)  = 0.0
180         VT20(1,3)  = 0.0
181         VT20(2,3)  = 1.0
182         VT20(3,3)  = 0.0
183         VT20(1,4)  = 0.0
184         VT20(2,4)  =-1.0
185         VT20(3,4)  = 0.0
186         VT20(1,5)  = 0.0
187         VT20(2,5)  = 0.0
188         VT20(3,5)  = 1.0
189         VT20(1,6)  = 0.0
190         VT20(2,6)  = 0.0
191         VT20(3,6)  =-1.0
192      END IF
193!
194      IF(NTSATM.EQ.8) THEN
195         VT20(1,1)  = 1.0
196         VT20(2,1)  = 1.0
197         VT20(3,1)  = 1.0
198         VT20(1,2)  =-1.0
199         VT20(2,2)  = 1.0
200         VT20(3,2)  = 1.0
201         VT20(1,3)  = 1.0
202         VT20(2,3)  =-1.0
203         VT20(3,3)  = 1.0
204         VT20(1,4)  = 1.0
205         VT20(2,4)  = 1.0
206         VT20(3,4)  =-1.0
207         VT20(1,5)  =-1.0
208         VT20(2,5)  =-1.0
209         VT20(3,5)  = 1.0
210         VT20(1,6)  =-1.0
211         VT20(2,6)  = 1.0
212         VT20(3,6)  =-1.0
213         VT20(1,7)  = 1.0
214         VT20(2,7)  =-1.0
215         VT20(3,7)  =-1.0
216         VT20(1,8)  =-1.0
217         VT20(2,8)  =-1.0
218         VT20(3,8)  =-1.0
219         CALL DSCAL(24,SQRT13,VT20,1)
220      END IF
221!
222      IF(NTSATM.EQ.12) THEN
223         VT20(1,1)  = 0.0
224         VT20(2,1)  = 1.0
225         VT20(3,1)  =    GOLD
226         VT20(1,2)  = 0.0
227         VT20(2,2)  = 1.0
228         VT20(3,2)  =   -GOLD
229         VT20(1,3)  = 0.0
230         VT20(2,3)  =-1.0
231         VT20(3,3)  =    GOLD
232         VT20(1,4)  = 0.0
233         VT20(2,4)  =-1.0
234         VT20(3,4)  =   -GOLD
235         VT20(3,5)  = 0.0
236         VT20(1,5)  = 1.0
237         VT20(2,5)  =    GOLD
238         VT20(3,6)  = 0.0
239         VT20(1,6)  = 1.0
240         VT20(2,6)  =   -GOLD
241         VT20(3,7)  = 0.0
242         VT20(1,7)  =-1.0
243         VT20(2,7)  =    GOLD
244         VT20(3,8)  = 0.0
245         VT20(1,8)  =-1.0
246         VT20(2,8)  =   -GOLD
247         VT20(2,9)  = 0.0
248         VT20(3,9)  = 1.0
249         VT20(1,9)  =    GOLD
250         VT20(2,10) = 0.0
251         VT20(3,10) = 1.0
252         VT20(1,10) =   -GOLD
253         VT20(2,11) = 0.0
254         VT20(3,11) =-1.0
255         VT20(1,11) =    GOLD
256         VT20(2,12) = 0.0
257         VT20(3,12) =-1.0
258         VT20(1,12) =   -GOLD
259         CALL DSCAL(36,0.525731112119134,VT20,1)
260      END IF
261!
262      IF(NTSATM.EQ.20) THEN
263         VT20(1,1)  = 1.0
264         VT20(2,1)  = 1.0
265         VT20(3,1)  = 1.0
266         VT20(1,2)  = 1.0
267         VT20(2,2)  = 1.0
268         VT20(3,2)  =-1.0
269         VT20(1,3)  = 1.0
270         VT20(2,3)  =-1.0
271         VT20(3,3)  = 1.0
272         VT20(1,4)  =-1.0
273         VT20(2,4)  = 1.0
274         VT20(3,4)  = 1.0
275         VT20(1,5)  = 1.0
276         VT20(2,5)  =-1.0
277         VT20(3,5)  =-1.0
278         VT20(1,6)  =-1.0
279         VT20(2,6)  = 1.0
280         VT20(3,6)  =-1.0
281         VT20(1,7)  =-1.0
282         VT20(2,7)  =-1.0
283         VT20(3,7)  = 1.0
284         VT20(1,8)  =-1.0
285         VT20(2,8)  =-1.0
286         VT20(3,8)  =-1.0
287         VT20(1,9)  = 0.0
288         VT20(2,9)  = ONEGOLD
289         VT20(3,9)  =    GOLD
290         VT20(1,10) = 0.0
291         VT20(2,10) = ONEGOLD
292         VT20(3,10) =   -GOLD
293         VT20(1,11) = 0.0
294         VT20(2,11) =-ONEGOLD
295         VT20(3,11) =    GOLD
296         VT20(1,12) = 0.0
297         VT20(2,12) =-ONEGOLD
298         VT20(3,12) =   -GOLD
299         VT20(3,13) = 0.0
300         VT20(1,13) = ONEGOLD
301         VT20(2,13) =    GOLD
302         VT20(3,14) = 0.0
303         VT20(1,14) = ONEGOLD
304         VT20(2,14) =   -GOLD
305         VT20(3,15) = 0.0
306         VT20(1,15) =-ONEGOLD
307         VT20(2,15) =    GOLD
308         VT20(3,16) = 0.0
309         VT20(1,16) =-ONEGOLD
310         VT20(2,16) =   -GOLD
311         VT20(2,17) = 0.0
312         VT20(3,17) = ONEGOLD
313         VT20(1,17) =    GOLD
314         VT20(2,18) = 0.0
315         VT20(3,18) = ONEGOLD
316         VT20(1,18) =   -GOLD
317         VT20(2,19) = 0.0
318         VT20(3,19) =-ONEGOLD
319         VT20(1,19) =    GOLD
320         VT20(2,20) = 0.0
321         VT20(3,20) =-ONEGOLD
322         VT20(1,20) =   -GOLD
323         CALL DSCAL(60,SQRT13,VT20,1)
324      END IF
325!
326      DO I = 1, NCENTS
327         INUC = INT(all_charges(I))
328         RFIX(I) = 2.400*aa2bohr
329
330         IF(INUC.EQ. 0) RFIX(I) = 0.001*aa2bohr
331         IF(INUC.EQ. 1) RFIX(I) = 1.400*aa2bohr
332         IF(INUC.EQ. 3) RFIX(I) = 1.400*aa2bohr
333         IF(INUC.EQ. 4) RFIX(I) = 1.400*aa2bohr
334         IF(INUC.EQ. 5) RFIX(I) = 1.400*aa2bohr
335         IF(INUC.EQ. 6) RFIX(I) = 2.100*aa2bohr
336         IF(INUC.EQ. 7) RFIX(I) = 2.000*aa2bohr
337         IF(INUC.EQ. 8) RFIX(I) = 1.900*aa2bohr
338         IF(INUC.EQ. 9) RFIX(I) = 1.800*aa2bohr
339         IF(INUC.EQ.10) RFIX(I) = 1.800*aa2bohr
340         IF(INUC.EQ.11) RFIX(I) = 1.800*aa2bohr
341         IF(INUC.EQ.12) RFIX(I) = 1.800*aa2bohr
342         IF(INUC.EQ.13) RFIX(I) = 1.800*aa2bohr
343         IF(INUC.EQ.14) RFIX(I) = 2.000*aa2bohr
344         IF(INUC.EQ.15) RFIX(I) = 2.200*aa2bohr
345         IF(INUC.EQ.16) RFIX(I) = 2.400*aa2bohr
346         IF(INUC.EQ.17) RFIX(I) = 2.760*aa2bohr
347         IF(INUC.EQ.18) RFIX(I) = 3.000*aa2bohr
348      END DO
349!
350!     -- NOTE: 'ME' STARTS AT 0 --
351!
352      xts = 0.0
353      yts = 0.0
354      zts = 0.0
355      afix = 0.0
356      idatom = 0
357      dai = 0.0
358      iddai = 0
359!      NFFTS   = ME*(MXFFTS/NPROC - 1)
360!      IPCOUNT = ME - 1
361      DO IFFAT = 1, NCENTS
362         IF(RFIX(IFFAT).LE.0.1) cycle
363!         IF(GOPARR) THEN
364!           IPCOUNT = IPCOUNT + 1
365!           IF(MOD(IPCOUNT,NPROC).NE.0) GOTO 500
366!         END IF
367         CALL FIXPVA2(IFFAT,all_centers,JVT1,CV,CDTST,AST,  &
368                    & XTS,YTS,ZTS,AFIX,RFIX,IDATOM,  &
369                    & DAI,IDDAI,DAIT,IDTMP,TMPTS,VT20)
370      end do
371!     Parallel code not available in dalton yet
372!     IF(GOPARR) THEN
373!        CALL DDI_GSUMF(2418,XTS   ,      MXFFTS)
374!        CALL DDI_GSUMF(2419,YTS   ,      MXFFTS)
375!        CALL DDI_GSUMF(2420,ZTS   ,      MXFFTS)
376!        CALL DDI_GSUMF(2421,AFIX  ,      MXFFTS)
377!        CALL DDI_GSUMI(2422,IDATOM,      MXFFTS)
378!        CALL DDI_GSUMF(2423,DAI   , 3*40*MXFFTS)
379!        CALL DDI_GSUMI(2424,IDDAI ,   41*MXFFTS)
380!     END IF
381!
382
383      KFFTS = 0
384      DO ITS=1,MXFFTS
385         IF(AFIX(ITS).LT.1.0e-4) cycle  ! 1.0e-04
386         KFFTS = KFFTS + 1
387         IF(KFFTS.GT.MXFFTS) THEN
388            ERROR STOP 'FIXTES: PLEASE INCREASE MXFFTS.'
389         END IF
390         XTS(KFFTS)      = XTS(ITS)
391         YTS(KFFTS)      = YTS(ITS)
392         ZTS(KFFTS)      = ZTS(ITS)
393         AFIX(KFFTS)     = AFIX(ITS)
394         IDATOM(KFFTS)   = IDATOM(ITS)
395         IDDAI(41,KFFTS) = IDDAI(41,ITS)
396         DO JJJ = 1, IDDAI(41,ITS)
397            IDDAI(JJJ,KFFTS) = IDDAI(JJJ,ITS)
398            DAI(1,JJJ,KFFTS) = DAI(1,JJJ,ITS)
399            DAI(2,JJJ,KFFTS) = DAI(2,JJJ,ITS)
400            DAI(3,JJJ,KFFTS) = DAI(3,JJJ,ITS)
401         END DO
402      end do
403      NFFTS = KFFTS   ! NFFTS IS GLOBAL
404
405!      FIXA = 0.0
406!      DO IFFTS=1,NFFTS
407!         FIXA = FIXA + AFIX(IFFTS)
408!      END DO
409!      write(luout,*) 'Surface point coordinates'
410!      do iffts= 1, nffts
411!         write(luout,*) XTS(IFFTS), YTS(IFFTS), ZTS(IFFTS)
412!      end do
413!      FIXA = FIXA*bohr2aa**2
414!
415!      write(luout,*) 'Surface area FIXA', FIXA, '/A**2'
416!      write(luout,*) 'Number of surface charges', NFFTS
417!      call openfile(surfile, surf, 'new', 'formatted')
418!      rewind(lu)
419!      write(surf,'(i6)') nffts
420!      write(surf,'(a)') 'AA'
421!      do i = 1, nffts
422!          write(surf,'(4f12.6)') XTS(I), YTS(I), ZTS(I), AFIX(I)
423!      end do
424      nsurp = nffts
425      allocate(Rsp(3,nsurp))
426      allocate(Sa(nsurp))
427      allocate(idatm(nsurp))
428      allocate(Radsp(nsurp))
429      do i = 1, nsurp
430         Sa(i) = AFIX(i)
431         Rsp(1,i) = XTS(i)
432         Rsp(2,i) = YTS(i)
433         Rsp(3,i) = ZTS(i)
434         idatm(i) = idatom(i)
435         Radsp(i) = rfix(idatom(i))
436      end do
437!      close(lu)
438
439end subroutine fixtes
440
441!-----------------------------------------------------------------------------
442
443subroutine fixpva2(iffat,cord,jvt1,cv,cdtst,ast,&
444                & xts,yts,zts,afix,rfix,idatom,&
445                & dai, iddai, dait,idtmp,tmpts,vt20)
446
447    use pelib_constants, only: pi
448
449    real(rp), dimension(:), intent(out) :: XTS, YTS, ZTS, AFIX, RFIX
450
451    real(rp), dimension(3,NCENTS) :: CORD
452    real(rp), dimension(960) :: AST
453    real(rp), dimension(3,960) :: CDTST
454    real(rp), dimension(3,40,960) :: DAIT
455    integer(ip), dimension(41,960) :: IDTMP
456    integer(ip), dimension(41) :: IDTMPTS
457    real(rp), dimension(6,60) :: JVT1
458    real(rp), dimension(122,3) :: CV
459    real(rp), dimension(3,20) :: VT20
460    real(rp), dimension(3,MXFFTS) :: TMPTS
461    real(rp) :: XII, YII, ZII, RII, AREA0, DUMMY
462    real(rp) :: P1X, P1Y, P1Z, P2X, P2Y, P2Z, P3X, P3Y, P3Z, P4X, P4Y, P4Z
463    real(rp) :: P5X, P5Y, P5Z, P6X, P6Y, P6Z, P7X, P7Y, P7Z
464    real(rp) :: PTS11, PTS21, PTS31, PTS12, PTS22, PTS32, PTS13, PTS23, PTS33
465    real(rp) :: DNORM4, DNORM5, DNORM6, DNORM7, FACTOR, SCALE4, SCALE5, SCALE6, SCALE7
466    integer(ip) :: IFFAT, KFFAT, ITS, KKK, MFFAT, III
467    integer(ip) :: N1, N2, N3
468    integer(ip) :: JJJ, JTS, KTS, LTS
469    real(rp) :: FOURPI = 4.0 * pi
470    real(rp) :: ONE3RD=1.0/3.0
471    real(rp), dimension(3,40,MXFFTS) :: DAI
472    integer(ip), dimension(MXFFTS) :: IDATOM
473    integer(ip), dimension(41,MXFFTS) :: IDDAI
474
475!     PARTITION OF THE CAVITY SURFACE INTO TESSERAE
476!     HUI LI, NOV 26, 2011
477
478     XII    = CORD(1,IFFAT)
479     YII    = CORD(2,IFFAT)
480     ZII    = CORD(3,IFFAT)
481     RII    = RFIX(IFFAT)
482     AREA0  = FOURPI*RII*RII/NTSATM
483     cdtst  = 0.0
484     ast    = 0.0
485     dait   = 0.0
486     idtmp  = 0
487
488
489
490!     --- USE 20 TESSERAE FOR EACH SPHERE ---
491     IF(NTSATM.EQ.20) THEN
492     DO ITS = 1, NTSATM
493
494        CDTST(1,ITS) = XII + VT20(1,ITS)*RII
495        CDTST(2,ITS) = YII + VT20(2,ITS)*RII
496        CDTST(3,ITS) = ZII + VT20(3,ITS)*RII
497
498        tmpts = 0.0
499        idtmpts = 0
500        KFFAT=IFFAT
501        AST(ITS) = AREA0
502        CALL FIXPVASWF(KFFAT,CORD,CDTST,ITS,AST(ITS),RFIX,&
503                     & TMPTS,IDTMPTS)
504        IF(AST(ITS).GT.0.9e-4) THEN
505           DO KKK = 1, IDTMPTS(41)
506              MFFAT = IDTMPTS(KKK)
507              DUMMY=ABS(TMPTS(1,MFFAT))+&
508                  & ABS(TMPTS(2,MFFAT))+&
509                  & ABS(TMPTS(3,MFFAT))
510              IF(DUMMY.GT.1.0e-5) THEN
511                 IDTMP(41,ITS) = IDTMP(41,ITS) + 1
512                 III = IDTMP(41,ITS)
513                 IDTMP(III,ITS)  = MFFAT
514                 DAIT(1,III,ITS) = TMPTS(1,MFFAT)
515                 DAIT(2,III,ITS) = TMPTS(2,MFFAT)
516                 DAIT(3,III,ITS) = TMPTS(3,MFFAT)
517              END IF
518           ENDDO
519        END IF
520     END DO
521     END IF
522
523
524
525!     --- USE 60 TESSERAE FOR EACH SPHERE ---
526     IF(NTSATM.EQ.60) THEN
527     DO ITS = 1, NTSATM
528
529        N1 = JVT1(1,ITS)
530        N2 = JVT1(2,ITS)
531        N3 = JVT1(3,ITS)
532        P4X = (CV(N1,1)+CV(N2,1)+CV(N3,1))*ONE3RD
533        P4Y = (CV(N1,2)+CV(N2,2)+CV(N3,2))*ONE3RD
534        P4Z = (CV(N1,3)+CV(N2,3)+CV(N3,3))*ONE3RD
535        DNORM4 = P4X**2+P4Y**2+P4Z**2
536        SCALE4 = 1.0/SQRT(DNORM4)
537        CDTST(1,ITS) = XII + P4X*SCALE4*RII
538        CDTST(2,ITS) = YII + P4Y*SCALE4*RII
539        CDTST(3,ITS) = ZII + P4Z*SCALE4*RII
540
541        tmpts = 0.0
542        idtmpts = 0
543        KFFAT=IFFAT
544        AST(ITS) = AREA0
545        CALL FIXPVASWF(KFFAT,CORD,CDTST,ITS,AST(ITS),RFIX,&
546                     & TMPTS,IDTMPTS)
547        IF(AST(ITS).GT.0.9D-04) THEN
548           DO KKK = 1, IDTMPTS(41)
549              MFFAT = IDTMPTS(KKK)
550              DUMMY=ABS(TMPTS(1,MFFAT))+&
551                  & ABS(TMPTS(2,MFFAT))+&
552                  & ABS(TMPTS(3,MFFAT))
553              IF(DUMMY.GT.1.0e-05) THEN
554                 IDTMP(41,ITS) = IDTMP(41,ITS) + 1
555                 III = IDTMP(41,ITS)
556                 IDTMP(III,ITS)  = MFFAT
557                 DAIT(1,III,ITS) = TMPTS(1,MFFAT)
558                 DAIT(2,III,ITS) = TMPTS(2,MFFAT)
559                 DAIT(3,III,ITS) = TMPTS(3,MFFAT)
560              END IF
561           ENDDO
562        END IF
563     end do
564     END IF
565
566
567
568!     --- USE 240 TESSERAE FOR EACH SPHERE ---
569     IF(NTSATM.EQ.240) THEN
570     DO KTS = 1, 60
571
572     DO JTS = 1, 4
573        ITS = (KTS-1)*4 + JTS
574        IF(JTS.EQ.1) THEN
575           N1 = JVT1(1,KTS)
576           N2 = JVT1(4,KTS)
577           N3 = JVT1(5,KTS)
578        END IF
579        IF(JTS.EQ.2) THEN
580           N1 = JVT1(6,KTS)
581           N2 = JVT1(4,KTS)
582           N3 = JVT1(5,KTS)
583        END IF
584        IF(JTS.EQ.3) THEN
585           N1 = JVT1(3,KTS)
586           N2 = JVT1(4,KTS)
587           N3 = JVT1(6,KTS)
588        END IF
589        IF(JTS.EQ.4) THEN
590           N1 = JVT1(2,KTS)
591           N2 = JVT1(6,KTS)
592           N3 = JVT1(5,KTS)
593        END IF
594        P4X = (CV(N1,1)+CV(N2,1)+CV(N3,1))*ONE3RD
595        P4Y = (CV(N1,2)+CV(N2,2)+CV(N3,2))*ONE3RD
596        P4Z = (CV(N1,3)+CV(N2,3)+CV(N3,3))*ONE3RD
597        DNORM4 = P4X**2+P4Y**2+P4Z**2
598        SCALE4 = 1.0/SQRT(DNORM4)
599        CDTST(1,ITS) = XII + P4X*SCALE4*RII
600        CDTST(2,ITS) = YII + P4Y*SCALE4*RII
601        CDTST(3,ITS) = ZII + P4Z*SCALE4*RII
602        KKK = MOD(ITS,4)
603        IF(KKK.EQ.1) FACTOR = 0.97527227808
604        IF(KKK.EQ.2) FACTOR = 1.04680294088
605        IF(KKK.EQ.3) FACTOR = 0.98896281888
606        IF(KKK.EQ.0) FACTOR = 0.98896281888
607
608        tmpts = 0.0
609        idtmpts = 0
610        KFFAT=IFFAT
611        AST(ITS) = AREA0*FACTOR
612        CALL FIXPVASWF(KFFAT,CORD,CDTST,ITS,AST(ITS),RFIX,&
613                     & TMPTS,IDTMPTS)
614        IF(AST(ITS).GT.0.9D-04) THEN
615           DO KKK = 1, IDTMPTS(41)
616              MFFAT = IDTMPTS(KKK)
617              DUMMY=ABS(TMPTS(1,MFFAT))+&
618                  & ABS(TMPTS(2,MFFAT))+&
619                  & ABS(TMPTS(3,MFFAT))
620              IF(DUMMY.GT.1.0e-05) THEN
621                 IDTMP(41,ITS) = IDTMP(41,ITS) + 1
622                 III = IDTMP(41,ITS)
623                 IDTMP(III,ITS)  = MFFAT
624                 DAIT(1,III,ITS) = TMPTS(1,MFFAT)
625                 DAIT(2,III,ITS) = TMPTS(2,MFFAT)
626                 DAIT(3,III,ITS) = TMPTS(3,MFFAT)
627              END IF
628           ENDDO
629        END IF
630     END DO
631     END DO
632     END IF
633
634
635
636!     --- USE 960 TESSERAE FOR EACH SPHERE ---
637     IF(NTSATM.EQ.960) THEN
638     DO KTS = 1, 60
639
640     DO JTS = 1, 4
641        IF(JTS.EQ.1) THEN
642           N1 = JVT1(1,KTS)
643           N2 = JVT1(4,KTS)
644           N3 = JVT1(5,KTS)
645        END IF
646        IF(JTS.EQ.2) THEN
647           N1 = JVT1(6,KTS)
648           N2 = JVT1(4,KTS)
649           N3 = JVT1(5,KTS)
650        END IF
651        IF(JTS.EQ.3) THEN
652           N1 = JVT1(3,KTS)
653           N2 = JVT1(4,KTS)
654           N3 = JVT1(6,KTS)
655        END IF
656        IF(JTS.EQ.4) THEN
657           N1 = JVT1(2,KTS)
658           N2 = JVT1(6,KTS)
659           N3 = JVT1(5,KTS)
660        END IF
661        P1X = CV(N1,1)
662        P1Y = CV(N1,2)
663        P1Z = CV(N1,3)
664        P2X = CV(N2,1)
665        P2Y = CV(N2,2)
666        P2Z = CV(N2,3)
667        P3X = CV(N3,1)
668        P3Y = CV(N3,2)
669        P3Z = CV(N3,3)
670!          COMPUTE THE COORDINATES OF POINTS 4, 5, 6
671!
672!                         1
673!
674!                      4     5
675!
676!                   3     6     2
677
678
679        P4X = (P1X+P3X)/2.0
680        P4Y = (P1Y+P3Y)/2.0
681        P4Z = (P1Z+P3Z)/2.0
682        P5X = (P1X+P2X)/2.0
683        P5Y = (P1Y+P2Y)/2.0
684        P5Z = (P1Z+P2Z)/2.0
685        P6X = (P2X+P3X)/2.0
686        P6Y = (P2Y+P3Y)/2.0
687        P6Z = (P2Z+P3Z)/2.0
688
689        DNORM4 = P4X**2+P4Y**2+P4Z**2
690        DNORM5 = P5X**2+P5Y**2+P5Z**2
691        DNORM6 = P6X**2+P6Y**2+P6Z**2
692        SCALE4 = 1.0/SQRT(DNORM4)
693        SCALE5 = 1.0/SQRT(DNORM5)
694        SCALE6 = 1.0/SQRT(DNORM6)
695        P4X = P4X*SCALE4
696        P4Y = P4Y*SCALE4
697        P4Z = P4Z*SCALE4
698        P5X = P5X*SCALE5
699        P5Y = P5Y*SCALE5
700        P5Z = P5Z*SCALE5
701        P6X = P6X*SCALE6
702        P6Y = P6Y*SCALE6
703        P6Z = P6Z*SCALE6
704
705    DO  LTS = 1, 4
706        ITS = ((KTS-1)*4 + JTS-1)*4 + LTS
707        IF(LTS.EQ.1) THEN
708          PTS11=P1X
709          PTS21=P1Y
710          PTS31=P1Z
711          PTS12=P4X
712          PTS22=P4Y
713          PTS32=P4Z
714          PTS13=P5X
715          PTS23=P5Y
716          PTS33=P5Z
717        ELSE IF(LTS.EQ.2) THEN
718          PTS11=P6X
719          PTS21=P6Y
720          PTS31=P6Z
721          PTS12=P4X
722          PTS22=P4Y
723          PTS32=P4Z
724          PTS13=P5X
725          PTS23=P5Y
726          PTS33=P5Z
727        ELSE IF(LTS.EQ.3) THEN
728          PTS11=P3X
729          PTS21=P3Y
730          PTS31=P3Z
731          PTS12=P4X
732          PTS22=P4Y
733          PTS32=P4Z
734          PTS13=P6X
735          PTS23=P6Y
736          PTS33=P6Z
737        ELSE IF(LTS.EQ.4) THEN
738          PTS11=P2X
739          PTS21=P2Y
740          PTS31=P2Z
741          PTS12=P6X
742          PTS22=P6Y
743          PTS32=P6Z
744          PTS13=P5X
745          PTS23=P5Y
746          PTS33=P5Z
747        END IF
748
749        P7X = (PTS11+PTS12+PTS13)*ONE3RD
750        P7Y = (PTS21+PTS22+PTS23)*ONE3RD
751        P7Z = (PTS31+PTS32+PTS33)*ONE3RD
752        DNORM7 = P7X**2+P7Y**2+P7Z**2
753        SCALE7 = 1.0/SQRT(DNORM7)
754        CDTST(1,ITS) = XII + P7X*SCALE7*RII
755        CDTST(2,ITS) = YII + P7Y*SCALE7*RII
756        CDTST(3,ITS) = ZII + P7Z*SCALE7*RII
757        KKK = MOD(ITS,16)
758        IF(KKK.EQ. 1) FACTOR = 0.96853843384
759        IF(KKK.EQ. 2) FACTOR = 0.98634758299
760        IF(KKK.EQ. 3) FACTOR = 0.97310155328
761        IF(KKK.EQ. 4) FACTOR = 0.97310155328
762        IF(KKK.EQ. 5) FACTOR = 1.04026074020
763        IF(KKK.EQ. 6) FACTOR = 1.05941468448
764        IF(KKK.EQ. 7) FACTOR = 1.04376817374
765        IF(KKK.EQ. 8) FACTOR = 1.04376817369
766        IF(KKK.EQ. 9) FACTOR = 0.98544774054
767        IF(KKK.EQ.10) FACTOR = 1.00020007354
768        IF(KKK.EQ.11) FACTOR = 0.98676349755
769        IF(KKK.EQ.12) FACTOR = 0.98343824680
770        IF(KKK.EQ.13) FACTOR = 0.98544774307
771        IF(KKK.EQ.14) FACTOR = 1.00020007615
772        IF(KKK.EQ.15) FACTOR = 0.98343824928
773        IF(KKK.EQ. 0) FACTOR = 0.98676350013
774
775        tmpts = 0.0
776        idtmpts = 0
777        KFFAT=IFFAT
778        AST(ITS) = AREA0*FACTOR
779        CALL FIXPVASWF(KFFAT,CORD,CDTST,ITS,AST(ITS),RFIX,&
780                     & TMPTS,IDTMPTS)
781       IF(AST(ITS).GT.0.9D-04) THEN
782           DO KKK = 1, IDTMPTS(41)
783              MFFAT = IDTMPTS(KKK)
784              DUMMY=ABS(TMPTS(1,MFFAT))+&
785                  & ABS(TMPTS(2,MFFAT))+&
786                  & ABS(TMPTS(3,MFFAT))
787              IF(DUMMY.GT.1.0e-05) THEN
788                 IDTMP(41,ITS) = IDTMP(41,ITS) + 1
789                 III = IDTMP(41,ITS)
790                 IDTMP(III,ITS)  = MFFAT
791                 DAIT(1,III,ITS) = TMPTS(1,MFFAT)
792                 DAIT(2,III,ITS) = TMPTS(2,MFFAT)
793                 DAIT(3,III,ITS) = TMPTS(3,MFFAT)
794              END IF
795           ENDDO
796        END IF
797     end do  ! lts = 1, 4
798     end do  ! jts = 1, 4
799     end do  ! kts = 1, 60
800     END IF
801
802
803     DO ITS=1,NTSATM
804        IF(AST(ITS).LT.0.95D-04) cycle      !  0.95D-04
805        NFFTS = NFFTS + 1
806        if (NFFTS .GT. MXFFTS) then
807           ERROR STOP 'FIXPVA2: PLEASE INCREASE MXFFTS.'
808        end if
809        XTS(NFFTS)      = CDTST(1,ITS)
810        YTS(NFFTS)      = CDTST(2,ITS)
811        ZTS(NFFTS)      = CDTST(3,ITS)
812        AFIX(NFFTS)     = AST(ITS)
813        IDATOM(NFFTS)   = IFFAT
814        IDDAI(41,NFFTS) = IDTMP(41,ITS)
815        DO JJJ = 1, IDTMP(41,ITS)
816           IDDAI(JJJ,NFFTS) = IDTMP(JJJ,ITS)
817           DAI(1,JJJ,NFFTS) = DAIT(1,JJJ,ITS)
818           DAI(2,JJJ,NFFTS) = DAIT(2,JJJ,ITS)
819           DAI(3,JJJ,NFFTS) = DAIT(3,JJJ,ITS)
820        ENDDO
821     END DO
822
823     RETURN
824
825end subroutine fixpva2
826
827!-----------------------------------------------------------------------------
828
829subroutine FIXPVASWF(IFFAT,CORD,CDTST,ITS,AREA,RFIX,TMP,IDTMPTS)
830
831    use pelib_constants
832
833   real(rp), intent(out) :: AREA
834   integer(ip), intent(in) :: iffat
835   real(rp), dimension(3,NCENTS) :: cord, tmp
836   real(rp), dimension(3,960) :: CDTST
837   real(rp), dimension(MXFFTS) :: RFIX
838   integer(ip), dimension(41) :: IDTMPTS
839   real(rp) :: ZERO=0.0
840   real(rp) :: ONE=1.0
841   real(rp) :: TWO=2.0
842   real(rp) :: PT5=0.5
843   real(rp) :: TEN=10.0
844   real(rp) :: DISM0, ONEDISM0, DISN1, DISN2, RA, RB
845   real(rp) :: X1, Y1, Z1, X2, Y2, Z2, X3, Y3, Z3, X5, Y5, Z5
846   real(rp) :: DISC2, ONEDISC2, DISC, ONEDISC, DISC3, DIS132, DIS13, DUM
847   real(rp) :: DISN, SWF2, DSWF2X, DSWF2Y, DSWF2Z
848   real(rp) :: DUW, DX5DX3, DY5DX3, DZ5DX3, DX5DY3, DY5DY3, DZ5DY3
849   real(rp) :: DX5DZ3, DY5DZ3, DZ5DZ3, DDISN2X, DDISN2Y, DDISN2Z
850   real(rp) :: SWF1, DSWF1X, DSWF1Y, DSWF1Z, DISD
851   real(rp) :: FB, FA, R1, R2, R3, R4, R5, R6, R7, R8, FAB, ONEFAB, DISM
852   real(rp) :: ONEDISM, ONEDISD, DDDIS
853   real(rp) :: CX3, CY3, CZ3, DX3, DY3, DZ3, R1X, R1Y, R1Z
854   real(rp) :: R2X, R2Y, R2Z, R3X, R3Y, R3Z, R4X, R4Y, R4Z, R5X, R5Y, R5Z
855   real(rp) :: R6X, R6Y, R6Z, R7X, R7Y, R7Z, R8X, R8Y, R8Z
856   real(rp) :: FABX, FABY, FABZ, TEMP, DDISM2X, DDISM2Y, DDISM2Z, DUMMY
857   integer(ip) :: JFFAT, KKK, ITS, KFFAT
858
859!
860!     DETERMINE THE AREA AND AREA DERIVATIVES FOR A TESSERA
861!     USING FIXPVA2
862!     HUI LI, NOV 27, 2011, LINCOLN
863!
864!     P1 = X1, Y1, Z1 (THE CURRENT TESSERA)
865!     P2 = X2, Y2, Z2 (THE CENTER OF THE SPHERE OF THE TESSERA)
866!     P3 = X3, Y3, Z3 (THE CENTER OF A NEIGHBOR SPHERE)
867!
868!
869      DISM0 = 0.30*aa2bohr
870      ONEDISM0 = ONE/DISM0
871      DISN1 = 0.70*aa2bohr
872      DISN2 = 1.20*aa2bohr    ! DISN2 MUST BE < RA
873!
874      X1 = CDTST(1,ITS)
875      Y1 = CDTST(2,ITS)
876      Z1 = CDTST(3,ITS)
877      X2 = CORD(1,IFFAT)
878      Y2 = CORD(2,IFFAT)
879      Z2 = CORD(3,IFFAT)
880      RA = RFIX(IFFAT)
881!
882      DO JFFAT=1,NCENTS
883!
884!        IF AREA=0, RETURN
885!        SKIP DISTANT SPHERES
886!
887         IF(JFFAT.EQ.IFFAT .OR. RFIX(JFFAT).LT.0.10) cycle
888         X3 = CORD(1,JFFAT)
889         Y3 = CORD(2,JFFAT)
890         Z3 = CORD(3,JFFAT)
891         RB = RFIX(JFFAT)
892         IF(ABS(X2-X3).GT.(RA+RB+DISN2)) cycle
893         IF(ABS(Y2-Y3).GT.(RA+RB+DISN2)) cycle
894         IF(ABS(Z2-Z3).GT.(RA+RB+DISN2)) cycle
895         DISC2 = (X2-X3)**2 + (Y2-Y3)**2 + (Z2-Z3)**2
896         ONEDISC2 = ONE/DISC2
897         IF(DISC2.GE.(RA+RB+DISN2)**2) cycle
898         IF(DISC2.LE.(RA-RB)**2 .AND. RB.LT.RA) cycle
899         IF(DISC2.LE.(RA-RB)**2 .AND. RB.GT.RA) THEN
900            AREA = ZERO
901            RETURN
902         END IF
903         DISC   = SQRT(DISC2)
904         ONEDISC= ONE/DISC
905         DISC3  = DISC*DISC2
906         DIS132 = (X1-X3)**2+(Y1-Y3)**2+(Z1-Z3)**2
907         DIS13  = SQRT(DIS132)
908!
909         DUM = RB*ONEDISC
910         X5  = X3 + (X2-X3)*DUM
911         Y5  = Y3 + (Y2-Y3)*DUM
912         Z5  = Z3 + (Z2-Z3)*DUM
913         DISN = SQRT((X1-X5)**2+(Y1-Y5)**2+(Z1-Z5)**2)
914         IF(DISN.GE.DISN2 .OR. DISC.LE.RB) THEN
915            SWF2   = ONE
916            DSWF2X = ZERO
917            DSWF2Y = ZERO
918            DSWF2Z = ZERO
919         ELSE IF (DISN.LE.DISN1) THEN
920            AREA   = ZERO
921            RETURN
922         ELSE
923            DUW  = (DISN**2 - DISN1**2)/(DISN2**2-DISN1**2)
924            SWF2 = 10.0*DUW**3 - 15.0*DUW**4 + 6.0*DUW**5
925            DUM    = RB/DISC3
926            DX5DX3 = ONE - RB*ONEDISC + DUM*(X2-X3)*(X2-X3)
927            DY5DX3 =                    DUM*(Y2-Y3)*(X2-X3)
928            DZ5DX3 =                    DUM*(Z2-Z3)*(X2-X3)
929            DX5DY3 =                    DUM*(X2-X3)*(Y2-Y3)
930            DY5DY3 = ONE - RB*ONEDISC + DUM*(Y2-Y3)*(Y2-Y3)
931            DZ5DY3 =                    DUM*(Z2-Z3)*(Y2-Y3)
932            DX5DZ3 =                    DUM*(X2-X3)*(Z2-Z3)
933            DY5DZ3 =                    DUM*(Y2-Y3)*(Z2-Z3)
934            DZ5DZ3 = ONE - RB*ONEDISC + DUM*(Z2-Z3)*(Z2-Z3)
935            DDISN2X=-TWO*((X1-X5)*DX5DX3+(Y1-Y5)*DY5DX3+(Z1-Z5)*DZ5DX3)
936            DDISN2Y=-TWO*((X1-X5)*DX5DY3+(Y1-Y5)*DY5DY3+(Z1-Z5)*DZ5DY3)
937            DDISN2Z=-TWO*((X1-X5)*DX5DZ3+(Y1-Y5)*DY5DZ3+(Z1-Z5)*DZ5DZ3)
938            DUM = (30.0*DUW**2-60.0*DUW**3+30.0*DUW**4)/&
939               &  (DISN2**2-DISN1**2)
940            DSWF2X = DUM*DDISN2X
941            DSWF2Y = DUM*DDISN2Y
942            DSWF2Z = DUM*DDISN2Z
943         END IF
944!
945!
946         IF(DISC.GE.(RA+RB)) THEN
947            SWF1   = ONE
948            DSWF1X = ZERO
949            DSWF1Y = ZERO
950            DSWF1Z = ZERO
951         ELSE
952            DISD = DIS13
953            FB =  RA**2 + DISC2 - DISD**2
954            FA =  RA**2 + DISC2 - RB**2
955            R1 =  RA    + DISC  + DISD
956            R2 =  RA    + DISC  - DISD
957            R3 =  RA    - DISC  + DISD
958            R4 = -RA    + DISC  + DISD
959            R5 =  RA    + DISC  + RB
960            R6 =  RA    + DISC  - RB
961            R7 =  RA    - DISC  + RB
962            R8 = -RA    + DISC  + RB
963            FAB = SQRT(ABS(R1*R2*R3*R4*R5*R6*R7*R8))
964            ONEFAB = ONE/FAB
965            DISM = SQRT(ABS(TWO*RA*RA - (FA*FB + FAB)*PT5*ONEDISC2))
966            IF(DISM.GE.DISM0) THEN
967               IF(DIS13.LT.RB) THEN
968                  AREA   = ZERO
969                  RETURN
970               END IF
971               SWF1   = ONE
972               DSWF1X = ZERO
973               DSWF1Y = ZERO
974               DSWF1Z = ZERO
975            ELSE
976               ONEDISM  = 1.0e+04
977               IF(DISM.GT.1.0e-04) ONEDISM  = ONE/DISM
978               ONEDISD  = ONE/DISD
979               DDDIS    = (DISM-DISM0)*ONEDISM0
980               SWF1     = PT5*DDDIS*DDDIS
981               IF(DIS13.GE.RB) SWF1 = ONE - SWF1
982               CX3 =  (X3-X2)*ONEDISC
983               CY3 =  (Y3-Y2)*ONEDISC
984               CZ3 =  (Z3-Z2)*ONEDISC
985               DX3 =  (X3-X1)*ONEDISD
986               DY3 =  (Y3-Y1)*ONEDISD
987               DZ3 =  (Z3-Z1)*ONEDISD
988               R1X =  CX3 + DX3
989               R1Y =  CY3 + DY3
990               R1Z =  CZ3 + DZ3
991               R2X =  CX3 - DX3
992               R2Y =  CY3 - DY3
993               R2Z =  CZ3 - DZ3
994               R3X = -CX3 + DX3
995               R3Y = -CY3 + DY3
996               R3Z = -CZ3 + DZ3
997               R4X =  CX3 + DX3
998               R4Y =  CY3 + DY3
999               R4Z =  CZ3 + DZ3
1000               R5X =  CX3
1001               R5Y =  CY3
1002               R5Z =  CZ3
1003               R6X =  CX3
1004               R6Y =  CY3
1005               R6Z =  CZ3
1006               R7X = -CX3
1007               R7Y = -CY3
1008               R7Z = -CZ3
1009               R8X =  CX3
1010               R8Y =  CY3
1011               R8Z =  CZ3
1012               FABX = (R1X*R2*R3*R4*R5*R6*R7*R8+&
1013                    &  R1*R2X*R3*R4*R5*R6*R7*R8+&
1014                    &  R1*R2*R3X*R4*R5*R6*R7*R8+&
1015                    &  R1*R2*R3*R4X*R5*R6*R7*R8+&
1016                    &  R1*R2*R3*R4*R5X*R6*R7*R8+&
1017                    &  R1*R2*R3*R4*R5*R6X*R7*R8+&
1018                    &  R1*R2*R3*R4*R5*R6*R7X*R8+&
1019                    &  R1*R2*R3*R4*R5*R6*R7*R8X)*ONEFAB*PT5
1020               FABY = (R1Y*R2*R3*R4*R5*R6*R7*R8+&
1021                    &  R1*R2Y*R3*R4*R5*R6*R7*R8+&
1022                    &  R1*R2*R3Y*R4*R5*R6*R7*R8+&
1023                    &  R1*R2*R3*R4Y*R5*R6*R7*R8+&
1024                    &  R1*R2*R3*R4*R5Y*R6*R7*R8+&
1025                    &  R1*R2*R3*R4*R5*R6Y*R7*R8+&
1026                    &  R1*R2*R3*R4*R5*R6*R7Y*R8+&
1027                    &  R1*R2*R3*R4*R5*R6*R7*R8Y)*ONEFAB*PT5
1028               FABZ = (R1Z*R2*R3*R4*R5*R6*R7*R8+&
1029                    &  R1*R2Z*R3*R4*R5*R6*R7*R8+&
1030                    &  R1*R2*R3Z*R4*R5*R6*R7*R8+&
1031                    &  R1*R2*R3*R4Z*R5*R6*R7*R8+&
1032                    &  R1*R2*R3*R4*R5Z*R6*R7*R8+&
1033                    &  R1*R2*R3*R4*R5*R6Z*R7*R8+&
1034                    &  R1*R2*R3*R4*R5*R6*R7Z*R8+&
1035                    &  R1*R2*R3*R4*R5*R6*R7*R8Z)*ONEFAB*PT5
1036               TEMP = (FA*FB + FAB)*ONEDISC2*ONEDISC2
1037               DDISM2X= (X3-X2)*TEMP&
1038                    &  -((X1-X2)*FA+(X3-X2)*FB)*ONEDISC2&
1039                    &  -FABX*PT5*ONEDISC2
1040               DDISM2Y= (Y3-Y2)*TEMP&
1041                    &  -((Y1-Y2)*FA+(Y3-Y2)*FB)*ONEDISC2&
1042                    &  -FABY*PT5*ONEDISC2
1043               DDISM2Z= (Z3-Z2)*TEMP&
1044                    &  -((Z1-Z2)*FA+(Z3-Z2)*FB)*ONEDISC2&
1045                    &  -FABZ*PT5*ONEDISC2
1046               DUM =    (DISM-DISM0)*ONEDISM0*ONEDISM0&
1047                    &   *PT5*ONEDISM
1048               IF(DIS13.GE.RB) DUM = -DUM
1049               DSWF1X = DUM*DDISM2X
1050               DSWF1Y = DUM*DDISM2Y
1051               DSWF1Z = DUM*DDISM2Z
1052               IF(DSWF1X.GT. TEN) DSWF1X =  TEN  ! AVOID SINGULARITY WHEN DISM=0
1053               IF(DSWF1X.LT.-TEN) DSWF1X = -TEN
1054               IF(DSWF1Y.GT. TEN) DSWF1Y =  TEN
1055               IF(DSWF1Y.LT.-TEN) DSWF1Y = -TEN
1056               IF(DSWF1Z.GT. TEN) DSWF1Z =  TEN
1057               IF(DSWF1Z.LT.-TEN) DSWF1Z = -TEN
1058            END IF
1059         END IF
1060!
1061!        - NOTE: LOOP 150 MEANS MULTI-SPHERE SCALING
1062!                AREA MUST BE SCALED AFTER AREA DERIVATIVES
1063!
1064         TMP(1,JFFAT) = AREA*(SWF1*DSWF2X + SWF2*DSWF1X)
1065         TMP(2,JFFAT) = AREA*(SWF1*DSWF2Y + SWF2*DSWF1Y)
1066         TMP(3,JFFAT) = AREA*(SWF1*DSWF2Z + SWF2*DSWF1Z)
1067         IF(JFFAT.GE.2) THEN
1068            DO KFFAT = 1, JFFAT-1
1069               TMP(1,KFFAT) = TMP(1,KFFAT)*SWF1*SWF2
1070               TMP(2,KFFAT) = TMP(2,KFFAT)*SWF1*SWF2
1071               TMP(3,KFFAT) = TMP(3,KFFAT)*SWF1*SWF2
1072            ENDDO
1073         END IF
1074         AREA = AREA*SWF1*SWF2
1075         IF(AREA.LT.0.9D-04) THEN   ! 0.9D-04
1076            AREA = ZERO
1077            RETURN
1078         END IF
1079         DUMMY=ABS(TMP(1,JFFAT))+ABS(TMP(2,JFFAT))+ABS(TMP(3,JFFAT))
1080         IF(DUMMY.LT.0.9D-05) cycle
1081         IDTMPTS(41)  = IDTMPTS(41) + 1
1082         KKK          = IDTMPTS(41)
1083         IF(KKK.EQ.39) THEN
1084            ERROR STOP 'FIXPVASWF: IDTMPTS EXCEEDED 40.'
1085         END IF
1086         IDTMPTS(KKK) = JFFAT
1087      end do ! jffat =1,NCENTS
1088!
1089!
1090      IDTMPTS(41)  = IDTMPTS(41) + 1
1091      KKK          = IDTMPTS(41)
1092      IDTMPTS(KKK) = IFFAT
1093      DO KFFAT = 1, NCENTS
1094         IF(KFFAT.NE.IFFAT) THEN
1095            TMP(1,IFFAT) = TMP(1,IFFAT) - TMP(1,KFFAT)
1096            TMP(2,IFFAT) = TMP(2,IFFAT) - TMP(2,KFFAT)
1097            TMP(3,IFFAT) = TMP(3,IFFAT) - TMP(3,KFFAT)
1098         END IF
1099      ENDDO
1100!
1101      RETURN
1102end subroutine FIXPVASWF
1103
1104end module pelib_cavity_generators
1105