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