1* 2*=============== DIPOLE =========================================== 3* 4* 9. Calculate dipole moment 5* -------------------------- 6 SUBROUTINE DIPOLE 7 include "prcm.h" 8* 9 FAK = 1.0/0.52917 10 FAC = 2.5418 11 N = 0 12 DO ITYP = 1,NTYPES 13 NSP = NSPEC(ITYP) 14 NSS = NSITS(ITYP) 15 FNSPI = 1.D0/DFLOAT(NSP) 16 DIPMOM = 0.D0 17 ISB = ISADR(ITYP)+1 18 ISE = ISADR(ITYP +1) 19 IBE = IADDR(ITYP)+1 20 IEN = IADDR(ITYP +1) 21 DO I = IBE,IEN 22 if(NSS.gt.1)then 23 DMX = 0.D0 24 DMY = 0.D0 25 DMZ = 0.D0 26 DO IS = ISB,ISE 27 N = N+1 28Calculate permanent dipole moment (magnitude) 29 DMX = DMX+CHARGE(IS)*WX(N) 30 DMY = DMY+CHARGE(IS)*WY(N) 31 DMZ = DMZ+CHARGE(IS)*WZ(N) 32 END DO! OF IS 33 DNORM = DSQRT(DMX**2+DMY**2+DMZ**2) 34 DIPMOM = DIPMOM+DNORM*FAC*FAK 35 IF(DNORM.GT.1.d-10) then 36 DPX(I) = DMX/DNORM 37 DPY(I) = DMY/DNORM 38 DPZ(I) = DMZ/DNORM 39 else 40 DPX(I)=0.d0 41 DPY(I)=0.d0 42 DPZ(I)=0.d0 43 end if 44 else 45 DIPMOM=0.d0 46 end if 47 END DO! OF I 48 DIPMOM = DIPMOM*FNSPI 49 IF(MOD(NSTEP,10).EQ.0.and.IPRINT.ge.8) THEN 50 IF(DABS(DIPMOM).GT.0.00001D0) THEN 51 PRINT "(' DIPOLE MOMENT(D):',F8.3,' FOR TYPE ',A6)" 52 X,DIPMOM,NAME(ITYP) 53 END IF 54 END IF 55 END DO! OF ITYP 56 RETURN 57 END 58* 59*=============== MPOLES ============================================= 60* 61 SUBROUTINE MPOLES(SITE,Q,DM,QM,OM,HM,NSB,NSE,NS) 62 IMPLICIT real*8 (A-H,O-Z) 63* 64 DIMENSION Q(NS),SITE(NS,3) 65 X,DM(3),QM(3,3),OM(3,3,3),HM(3,3,3,3) 66* 67 FAK = 1.0/0.52917 68 DO 10 IS = NSB,NSE 69 DO 10 J = 1,3 70 DM(J) = 0.D0 71 DO 10 K = 1,3 72 QM(J,K) = 0.D0 73 DO 10 L = 1,3 74 OM(J,K,L) = 0.D0 75 DO 10 M = 1,3 76 HM(J,K,L,M) = 0.D0 77 10 CONTINUE 78* 79 DO 20 IS = NSB,NSE 80 DO 20 J = 1,3 81 DM(J) = DM(J)+Q(IS)*SITE(IS,J)*FAK 82 DO 20 K = 1,3 83 QM(J,K) = QM(J,K)+Q(IS)*SITE(IS,J)*SITE(IS,K)*FAK**2 84 DO 20 L = 1,3 85 OM(J,K,L) = OM(J,K,L)+Q(IS)*SITE(IS,J)* 86 X SITE(IS,K)*SITE(IS,L)*FAK**3 87 DO 20 M = 1,3 88 HM(J,K,L,M) = HM(J,K,L,M)+Q(IS)*SITE(IS,J)* 89 X SITE(IS,K)*SITE(IS,L)*SITE(IS,M)*FAK**4 90 20 CONTINUE 91 RETURN 92 END 93 94* 95*=============== GETMOI ======================================== 96* 97 98 SUBROUTINE GETMOI(X,Y,Z,M,MOMENT,NSS) 99 IMPLICIT real*8 (A-H,O-Z) 100 real*8 M,MOMENT 101 DIMENSION X(NSS),Y(NSS),Z(NSS),M(NSS) 102 DIMENSION MOMENT(3,3),R(3) 103* 104 DO INERT = 1,3 105 DO IA = 1,3 106 TT = 0.D0 107 MOMENT(INERT,IA) = 0.D0 108 IF(INERT.EQ.IA) TT = 1.D0 109 DO I = 1,NSS 110 R(1) = X(I) 111 R(2) = Y(I) 112 R(3) = Z(I) 113 SQS = X(I)**2+Y(I)**2+Z(I)**2 114 MOMENT(INERT,IA) = MOMENT(INERT,IA)+M(I)*(TT*SQS-R(INERT)*R(IA)) 115 END DO! OF I 116 END DO! OF IA 117 END DO! OF INERT 118* 119 RETURN 120 END 121* 122*=============== GETROT ========================================== 123* 124 SUBROUTINE GETROT 125* 126 include"prcm.h" 127* 128 DIMENSION RX(NS),RY(NS),RZ(NS),FMASS(NS) 129 DIMENSION FMOI(3,3) 130 DIMENSION G(3),H(3) 131* 132 N = 0 133 DO ITYP = 1,NTYPES 134 NSS = NSITS (ITYP) 135 NSP = NSPEC (ITYP) 136 ISB = ISADR (ITYP)+1 137 ISE = ISADR (ITYP +1) 138 JBE = IADDR (ITYP)+1 139 JEN = IADDR (ITYP +1) 140 DO J = JBE,JEN 141 I = 0 142 DO IS = ISB,ISE 143 I = I+1 144 N = N+1 145 RX (I) = WX(N) 146 RY (I) = WY(N) 147 RZ (I) = WZ(N) 148 FMASS(I) = MASSD(IS) 149 END DO! OF IS 150* 151 CALL GETMOI(RX,RY,RZ,FMASS,FMOI,NSS) 152* 153 CALL HH3BY3(FMOI,G,H) 154* 155 DO K = 1,3 156 DO L = 1,3 157 RMX(J,K,L) = FMOI(K,L) 158 END DO! OF K 159 END DO! OF L 160* 161 END DO! OF J 162* 163 END DO! OF ITYP 164* 165 RETURN 166 END 167* 168*========== HH3BY3 ============================================ 169* 170 SUBROUTINE HH3BY3(A,D,E) 171* 172 IMPLICIT real*8 (A-H,O-Z) 173* 174 DIMENSION A(3,3),D(3),E(3) 175 data ICOUNT/0/ 176* 177 H = 0. 178 SCALE = ABS(A(3,1))+ABS(A(3,2)) 179 IF(SCALE.EQ.0) THEN 180 E(3) = A(3,2) 181 ELSE 182 A(3,1) = A(3,1)/SCALE 183 A(3,2) = A(3,2)/SCALE 184 H = A(3,1)**2+A(3,2)**2+H 185 F = A(3,2) 186 G =-SIGN(SQRT(H),F) 187 E(3) = SCALE*G 188 H = H-F*G 189 A(3,2) = F-G 190* 191 A(1,3) = A(3,1)/H 192 G = A(1,1)*A(3,1)+A(2,1)*A(3,2) 193 E(1) = G/H 194 F = E(1)*A(3,1) 195* 196 A(2,3) = A(3,2)/H 197 G = A(2,1)*A(3,1)+A(2,2)*A(3,2) 198 E(2) = G/H 199 F = F+E(2)*A(3,2) 200* 201 HH = F/(H+H) 202 F = A(3,1) 203 G = E(1)-HH*F 204 E(1) = G 205 A(1,1) = A(1,1)-F*E(1)-G*A(3,1) 206 F = A(3,2) 207 G = E(2)-HH*F 208 E(2) = G 209 A(2,1) = A(2,1)-F*E(1)-G*A(3,1) 210 A(2,2) = A(2,2)-F*E(2)-G*A(3,2) 211 END IF 212 D(3) = H 213* 214 SCALE = 0. 215 E(2) = A(2,1) 216 D(2) = 0. 217* 218 E(1) = 0. 219 D(1) = A(1,1) 220 A(1,1) = 1. 221 IF(D(2).NE.0) THEN 222 G = A(2,1)*A(1,1) 223 A(1,1) = A(1,1)-A(1,2)*G 224 END IF 225 D(2) = A(2,2) 226 A(2,2) = 1. 227 A(2,1) = 0. 228 A(1,2) = 0. 229 IF(D(3).NE.0) THEN 230 G = A(3,1)*A(1,1)+A(3,2)*A(2,1) 231 A(1,1) = A(1,1)-A(1,3)*G 232 A(2,1) = A(2,1)-A(2,3)*G 233 G = A(3,1)*A(1,2)+A(3,2)*A(2,2) 234 A(1,2) = A(1,2)-A(1,3)*G 235 A(2,2) = A(2,2)-A(2,3)*G 236 END IF 237 D(3) = A(3,3) 238 A(3,3) = 1. 239 A(3,1) = 0. 240 A(1,3) = 0. 241 A(3,2) = 0. 242 A(2,3) = 0. 243* 244*Diagonalize it: 245* 246 E(1) = E(2) 247 E(2) = E(3) 248 E(3) = 0. 249 DO 15 L = 1,3 250 ITER = 0 251 1 CONTINUE 252 DO 12 M = L,2 253 DD = ABS(D(M))+ABS(D(M+1)) 254 IF(ABS(E(M))+DD.EQ.DD) GO TO 2 255 12 CONTINUE 256 M = 3 257 2 CONTINUE 258 IF(M.NE.L) THEN 259 IF(ITER.EQ.30) then 260 write(*,*) '!!! TOO MANY ITERATIONS! in HH3BY3' 261 ICOUNT=ICOUNT+1 262 if(ICOUNT.gt.1000)stop 263 return 264 end if 265 ITER = ITER+1 266 G = (D(L+1)-D(L))/(2.*E(L)) 267 R = SQRT(G**2+1.) 268 G = D(M)-D(L)+E(L)/(G+SIGN(R,G)) 269 S = 1. 270 C = 1. 271 P = 0. 272 DO 14 I = M-1,L,-1 273 F = S*E(I) 274 B = C*E(I) 275 IF(ABS(F).GE.ABS(G)) THEN 276 C = G/F 277 R = SQRT(C**2+1.) 278 E(I+1) = F*R 279 S = 1./R 280 C = C*S 281 ELSE 282 S = F/G 283 R = SQRT(S**2+1.) 284 E(I+1) = G*R 285 C = 1./R 286 S = S*C 287 END IF 288 G = D(I+1)-P 289 R =(D(I)-G)*S+2.*C*B 290 P = S*R 291 D(I+1) = G+P 292 G = C*R-B 293 F = A(1,I+1) 294 A(1,I+1)= S*A(1,I)+C*F 295 A(1,I) = C*A(1,I)-S*F 296 F = A(2,I+1) 297 A(2,I+1)= S*A(2,I)+C*F 298 A(2,I) = C*A(2,I)-S*F 299 F = A(3,I+1) 300 A(3,I+1)= S*A(3,I)+C*F 301 A(3,I) = C*A(3,I)-S*F 302 14 CONTINUE 303 D(L) = D(L)-P 304 E(L) = G 305 E(M) = 0. 306 GO TO 1 307 END IF 308 15 CONTINUE 309 RETURN 310 END 311* 312*========= ROTMAT ================================================= 313* 314 SUBROUTINE ROTMAT(XX,YY,ZZ,IBG,N,NBTOP) 315* 316 include "prcm.h" 317* 318 DIMENSION XX(*),YY(*),ZZ(*) 319 LOGICAL NBTOP 320* 321 IBEG = IBG 322 IEND = IBG+N-1 323* 324 IF(NBTOP) THEN ! from box to principal: 325 DO I = IBEG,IEND 326 AQ = XX(I) 327 BQ = YY(I) 328 CQ = ZZ(I) 329 XX(I) = RMX(I,1,1)*AQ+RMX(I,2,1)*BQ+RMX(I,3,1)*CQ 330 YY(I) = RMX(I,1,2)*AQ+RMX(I,2,2)*BQ+RMX(I,3,2)*CQ 331 ZZ(I) = RMX(I,1,3)*AQ+RMX(I,2,3)*BQ+RMX(I,3,3)*CQ 332 END DO! OF I 333 ELSE ! from principal to box: 334 DO I = IBEG,IEND 335 AQ = XX(I) 336 BQ = YY(I) 337 CQ = ZZ(I) 338 XX(I) = RMX(I,1,1)*AQ+RMX(I,1,2)*BQ+RMX(I,1,3)*CQ 339 YY(I) = RMX(I,2,1)*AQ+RMX(I,2,2)*BQ+RMX(I,2,3)*CQ 340 ZZ(I) = RMX(I,3,1)*AQ+RMX(I,3,2)*BQ+RMX(I,3,3)*CQ 341 END DO! OF I 342 END IF 343* 344 RETURN 345 END 346* 347*============ FMELT ================================================ 348* 349 SUBROUTINE FMELT(X0,Y0,Z0,NSP) 350 IMPLICIT real*8 (A-H,O-Z) 351 DIMENSION X0(*),Y0(*),Z0(*) 352 DATA PI/3.1415926536D0/ 353 CN=PI*(2*NSP)**(1.D0/3.D0) 354 AB=0.D0 355 BC=0.D0 356 CD=0.D0 357 DE=0.D0 358 EF=0.D0 359 FG=0.D0 360 DO 10 I=1,NSP 361 AB=AB+DCOS(CN*X0(I)) 362 BC=BC+DSIN(CN*X0(I)) 363 CD=CD+DCOS(CN*Y0(I)) 364 DE=DE+DSIN(CN*Y0(I)) 365 EF=EF+DCOS(CN*Z0(I)) 366 FG=FG+DSIN(CN*Z0(I)) 367 10 CONTINUE 368 FMLT=(AB*AB+BC*BC+CD*CD+DE*DE+EF*EF+FG*FG)/DFLOAT(3*NSP) 369 PRINT "(/1X,' ***** MELTING FACTOR: ',F12.4/)",FMLT 370 RETURN 371 END 372* 373* Do read(STR,*)(NAME(I),I=1,N) where STR and NAME are characters 374* ------------------------------------------------------------------- 375C (only because some stupid compilers do not understand this) 376C========= PARSES ============================ 377C 378 subroutine PARSES(STR,CH,N,L,ie) 379 character STR(L),CHR 380 character*32 CH(N) 381 logical LI 382C LCH=32 is length of NAMOL in prcm.h 383 LCH=32 384 do I=1,N 385 do J=1,LCH 386 CH(I)(J:J)=' ' 387 end do 388 end do 389 LI=.false. 390 M=0 391 ie=0 392 do I=1,L 393 CHR=STR(I) 394 ICHR=ichar(CHR) 395 if(LI)then 396C Delimiters - space, null and tab 397 if(CHR.eq.' '.or.ICHR.eq.0.or.ICHR.eq.9)then 398 LI=.false. 399 else 400 IC=IC+1 401 if(IC.le.LCH)CH(M)(IC:IC)=CHR 402 end if 403 else 404 if(CHR.ne.' '.and.ICHR.ne.0.and.ICHR.ne.9)then 405 M=M+1 406 if(M.gt.N)return 407 CH(M)(1:1)=CHR 408 IC=1 409 LI=.true. 410 end if 411 end if 412 end do 413 if(M.lt.N)ie=1 414 return 415 end 416* 417*=============== NEXTITEM ================================= 418* 419 function NEXTITEM(STR) 420 character*128 STR 421 do I=1,128 422 if(STR(I:I).ne.' ')go to 10 423 end do 424 NEXTITEM=0 425 return 426 10 do J=I,128 427 if(STR(J:J).eq.' ')go to 20 428 end do 429 NEXTITEM=128 430 return 431 20 do I=J,128 432 if(STR(I:I).ne.' ')go to 30 433 end do 434 30 NEXTITEM=I 435 return 436 end 437* 438*======================== TRANSRDF ======================== 439* 440 subroutine TRANSRDF(NS,MXGRS,MAST,TASKID,NGRI,IGRI,MGRI,IAUX) 441 include "mpif.h" 442 integer NGRI(NS),IGRI(NS,MXGRS),MGRI(NS,MXGRS),IAUX(NS),TASKID 443 call MPI_BCAST(NGRI,NS,MPI_INTEGER,MAST,MPI_COMM_WORLD,ie) 444 do I=1,MXGRS 445 if(TASKID.eq.MAST)then 446 do J=1,NS 447 IAUX(J)=IGRI(J,I) 448 end do 449 end if 450 call MPI_BCAST(IAUX,NS,MPI_INTEGER,MAST,MPI_COMM_WORLD,ie) 451 if(TASKID.eq.MAST)then 452 do J=1,NS 453 IAUX(J)=MGRI(J,I) 454 end do 455 else 456 do J=1,NS 457 IGRI(J,I)=IAUX(J) 458 end do 459 end if 460 call MPI_BCAST(IAUX,NS,MPI_INTEGER,MAST,MPI_COMM_WORLD,ie) 461 if(TASKID.ne.MAST)then 462 do J=1,NS 463 MGRI(J,I)=IAUX(J) 464 end do 465 end if 466 end do 467 return 468 end 469 470