1 subroutine rotate(v,w,g,x,y) 2c 3c $Id$ 4c 5 implicit none 6c 7 real*8 v(3),w(3),x(3),y(3),xx(3),t(3,3) 8 real*8 small,pi,r,a,b,ca,cb,cg,sa,sb,sg,g 9 integer i 10 parameter (small=1.0d-24) 11c 12c rotation with angle g around vector from v to w 13c of point x giving result in y 14c 15 if(abs(w(2)-v(2)).lt.small) then 16 if(abs(w(1)-v(1)).lt.small) then 17 a=0.0d0 18 else 19 if(w(1)-v(1).gt.0.0d0) then 20 a=2.0d0*datan(1.0d0) 21 else 22 a=(-2.0d0)*datan(1.0d0) 23 endif 24 endif 25 else 26 a=atan(abs(w(1)-v(1))/abs(w(2)-v(2))) 27 pi=4.0d0*atan(1.0d0) 28 if(w(1)-v(1).gt.0.0d0.and.w(2)-v(2).lt.0.0d0) a=pi-a 29 if(w(1)-v(1).lt.0.0d0.and.w(2)-v(2).gt.0.0d0) a=-a 30 if(w(1)-v(1).lt.0.0d0.and.w(2)-v(2).lt.0.0d0) a=pi+a 31 endif 32 r=0.0d0 33 do 1 i=1,3 34 r=r+(w(i)-v(i))**2 35 xx(i)=x(i)-v(i) 36 1 continue 37 if(r.lt.small) then 38 y(1)=x(1) 39 y(2)=x(2) 40 y(3)=x(3) 41 return 42 endif 43 b=acos((w(3)-v(3))/sqrt(r)) 44 sa=sin(a) 45 ca=cos(a) 46 sb=sin(b) 47 cb=cos(b) 48 sg=sin(g) 49 cg=cos(g) 50 t(1,1)=ca*ca*cg-sa*ca*cb*sg+sa*ca*cb*sg 51 + +sa*sa*cb*cb*cg+sa*sa*sb*sb 52 t(1,2)=(-sa)*ca*cg-ca*ca*cb*sg-sa*sa*cb*sg 53 + +sa*ca*cb*cb*cg+sa*ca*sb*sb 54 t(1,3)=ca*sb*sg-sa*sb*cb*cg+sa*sb*cb 55 t(2,1)=(-sa)*ca*cg+sa*sa*cb*sg+ca*ca*cb*sg 56 + +sa*ca*cb*cb*cg+sa*ca*sb*sb 57 t(2,2)=sa*sa*cg+sa*ca*cb*sg-sa*ca*cb*sg 58 + +ca*ca*cb*cb*cg+ca*ca*sb*sb 59 t(2,3)=(-sa)*sb*sg-ca*sb*cb*cg+ca*sb*cb 60 t(3,1)=(-ca)*sb*sg-sa*sb*cb*cg+sa*sb*cb 61 t(3,2)=sa*sb*sg-ca*sb*cb*cg+ca*sb*cb 62 t(3,3)=sb*sb*cg+cb*cb 63 do 2 i=1,3 64 y(i)=xx(1)*t(i,1)+xx(2)*t(i,2)+xx(3)*t(i,3)+v(i) 65 2 continue 66 return 67 end 68 subroutine super(x,nx,mx,y,w,ws,sdev,ny,my,mod,rms0,rms1, 69 + xw,mw,nw,ma,na,lpara) 70c 71c superimpose x(1:n,1:3) onto y(1:n,1:3) 72c 73 implicit none 74c 75#include "msgids.fh" 76#include "global.fh" 77c 78 real*8 zero 79 parameter(zero=0.0d0) 80c 81 integer nx,mx,ny,my,mw,nw,ma,na 82 real*8 x(mx,3),y(my,3),w(my),ws(my),sdev(my),xw(mw,ma,3),rms0,rms1 83 logical mod,lpara 84c 85 integer i,j,k,l,nr 86 real*8 u(3,3),c(4,4),q(4),b(4),v(4,4),wnorm,wnorms 87 real*8 xt(3),yt(3),sd 88c 89c write(*,'(a,6f12.6)') 'Super in ',(x(1,i),i=1,3),(x(nx,i),i=1,3) 90 wnorm=zero 91 wnorms=zero 92 do 1 j=1,3 93 xt(j)=zero 94 yt(j)=zero 95 1 continue 96 do 2 k=1,nx 97 do 3 j=1,3 98 xt(j)=xt(j)+w(k)*x(k,j) 99 3 continue 100 2 continue 101 if(lpara) call ga_dgop(mag_d08,xt,3,'+') 102 do 4 k=1,ny 103 wnorm=wnorm+w(k) 104 wnorms=wnorms+ws(k) 105 do 5 j=1,3 106 yt(j)=yt(j)+w(k)*y(k,j) 107 5 continue 108 4 continue 109c 110 rms0=zero 111 do 6 j=1,3 112 xt(j)=xt(j)/wnorm 113 yt(j)=yt(j)/wnorm 114 do 7 k=1,nx 115 rms0=rms0+ws(k)*(x(k,j)-xt(j)-y(k,j)+yt(j))**2 116 7 continue 117 6 continue 118 if(lpara) call ga_dgop(mag_d09,rms0,1,'+') 119 rms0=sqrt(rms0/wnorms) 120c 121 do 8 i=1,3 122 do 9 j=1,3 123 u(i,j)=zero 124 do 10 k=1,nx 125 u(i,j)=u(i,j)+w(k)*(x(k,i)-xt(i))*(y(k,j)-yt(j)) 126 10 continue 127 9 continue 128 8 continue 129 if(lpara) call ga_dgop(mag_d10,u,9,'+') 130c 131 c(1,1)=u(1,1)+u(2,2)+u(3,3) 132 c(1,2)=u(3,2)-u(2,3) 133 c(1,3)=u(1,3)-u(3,1) 134 c(1,4)=u(2,1)-u(1,2) 135 c(2,2)=u(1,1)-u(2,2)-u(3,3) 136 c(2,3)=u(1,2)+u(2,1) 137 c(2,4)=u(3,1)+u(1,3) 138 c(3,3)=u(2,2)-u(3,3)-u(1,1) 139 c(3,4)=u(2,3)+u(3,2) 140 c(4,4)=u(3,3)-u(1,1)-u(2,2) 141c 142 do 11 j=1,3 143 do 12 i=j+1,4 144 c(i,j)=c(j,i) 145 12 continue 146 11 continue 147c 148 call md_jacobi(c,4,4,b,v,nr) 149c 150 do 13 i=1,4 151 q(i)=v(i,4) 152 13 continue 153c 154 u(1,1)=q(1)*q(1)+q(2)*q(2)-q(3)*q(3)-q(4)*q(4) 155 u(1,2)=2.0d0*(q(3)*q(2)+q(1)*q(4)) 156 u(1,3)=2.0d0*(q(4)*q(2)-q(1)*q(3)) 157 u(2,1)=2.0d0*(q(2)*q(3)-q(1)*q(4)) 158 u(2,2)=q(1)*q(1)-q(2)*q(2)+q(3)*q(3)-q(4)*q(4) 159 u(2,3)=2.0d0*(q(4)*q(3)+q(1)*q(2)) 160 u(3,1)=2.0d0*(q(2)*q(4)+q(1)*q(3)) 161 u(3,2)=2.0d0*(q(3)*q(4)-q(1)*q(2)) 162 u(3,3)=q(1)*q(1)-q(2)*q(2)-q(3)*q(3)+q(4)*q(4) 163c 164 rms1=zero 165 do 14 k=1,nx 166 do 15 i=1,3 167 q(i)=0.0d0 168 do 16 j=1,3 169 q(i)=q(i)+u(i,j)*(x(k,j)-xt(j)) 170 16 continue 171 15 continue 172 if(mod) then 173 do 17 i=1,3 174 x(k,i)=q(i)+yt(i) 175 17 continue 176 endif 177 sd=0.0d0 178 do 18 j=1,3 179 sd=sd+(q(j)-y(k,j)+yt(j))**2 180 18 continue 181 sdev(k)=sdev(k)+sd 182 rms1=rms1+ws(k)*sd 183 14 continue 184 if(lpara) call ga_dgop(mag_d11,rms1,1,'+') 185 rms1=sqrt(rms1/wnorms) 186c 187 if(mod.and.nw.gt.0) then 188 do 19 l=1,nw 189 do 20 k=1,na 190 do 21 i=1,3 191 q(i)=0.0d0 192 do 22 j=1,3 193 q(i)=q(i)+u(i,j)*(xw(l,k,j)-xt(j)) 194 22 continue 195 21 continue 196 do 23 i=1,3 197 xw(l,k,i)=q(i)+yt(i) 198 23 continue 199 20 continue 200 19 continue 201 endif 202c 203c write(*,'(a,6f12.6)') 'Super out',(x(1,i),i=1,3),(x(nx,i),i=1,3) 204 return 205 end 206 subroutine super2(x,ix,nx,mx,y,w,ws,sdev,ny,my,mod,rms0,rms1, 207 + xw,mw,nw,ma,na,lpara) 208c 209c superimpose x(1:n,1:3) onto y(1:n,1:3) 210c 211 implicit none 212c 213#include "msgids.fh" 214#include "global.fh" 215c 216 real*8 zero 217 parameter(zero=0.0d0) 218c 219 integer nx,mx,ny,my,mw,nw,ma,na 220 integer ix(mx) 221 real*8 x(mx,3),y(my,3),w(my),ws(my),sdev(my),xw(mw,ma,3),rms0,rms1 222 logical mod,lpara 223c 224 integer i,j,k,l,nr 225 real*8 u(3,3),c(4,4),q(4),b(4),v(4,4),wnorm,wnorms 226 real*8 xt(3),yt(3),sd 227c 228c write(*,'(a,6f12.6)') 'Super in ',(x(1,i),i=1,3),(x(nx,i),i=1,3) 229 wnorm=zero 230 wnorms=zero 231 do 1 j=1,3 232 xt(j)=zero 233 yt(j)=zero 234 1 continue 235 do 2 k=1,nx 236 do 3 j=1,3 237 xt(j)=xt(j)+w(ix(k))*x(k,j) 238 3 continue 239 2 continue 240 if(lpara) call ga_dgop(mag_d08,xt,3,'+') 241 do 4 k=1,ny 242 wnorm=wnorm+w(k) 243 wnorms=wnorms+ws(k) 244 do 5 j=1,3 245 yt(j)=yt(j)+w(k)*y(k,j) 246 5 continue 247 4 continue 248c 249 rms0=zero 250 do 6 j=1,3 251 xt(j)=xt(j)/wnorm 252 yt(j)=yt(j)/wnorm 253 do 7 k=1,nx 254 rms0=rms0+ws(ix(k))*(x(k,j)-xt(j)-y(ix(k),j)+yt(j))**2 255 7 continue 256 6 continue 257 if(lpara) call ga_dgop(mag_d09,rms0,1,'+') 258 rms0=sqrt(rms0/wnorms) 259c 260 do 8 i=1,3 261 do 9 j=1,3 262 u(i,j)=zero 263 do 10 k=1,nx 264 u(i,j)=u(i,j)+w(ix(k))*(x(k,i)-xt(i))*(y(ix(k),j)-yt(j)) 265 10 continue 266 9 continue 267 8 continue 268 if(lpara) call ga_dgop(mag_d10,u,9,'+') 269c 270 c(1,1)=u(1,1)+u(2,2)+u(3,3) 271 c(1,2)=u(3,2)-u(2,3) 272 c(1,3)=u(1,3)-u(3,1) 273 c(1,4)=u(2,1)-u(1,2) 274 c(2,2)=u(1,1)-u(2,2)-u(3,3) 275 c(2,3)=u(1,2)+u(2,1) 276 c(2,4)=u(3,1)+u(1,3) 277 c(3,3)=u(2,2)-u(3,3)-u(1,1) 278 c(3,4)=u(2,3)+u(3,2) 279 c(4,4)=u(3,3)-u(1,1)-u(2,2) 280c 281 do 11 j=1,3 282 do 12 i=j+1,4 283 c(i,j)=c(j,i) 284 12 continue 285 11 continue 286c 287 call md_jacobi(c,4,4,b,v,nr) 288c 289 do 13 i=1,4 290 q(i)=v(i,4) 291 13 continue 292c 293 u(1,1)=q(1)*q(1)+q(2)*q(2)-q(3)*q(3)-q(4)*q(4) 294 u(1,2)=2.0d0*(q(3)*q(2)+q(1)*q(4)) 295 u(1,3)=2.0d0*(q(4)*q(2)-q(1)*q(3)) 296 u(2,1)=2.0d0*(q(2)*q(3)-q(1)*q(4)) 297 u(2,2)=q(1)*q(1)-q(2)*q(2)+q(3)*q(3)-q(4)*q(4) 298 u(2,3)=2.0d0*(q(4)*q(3)+q(1)*q(2)) 299 u(3,1)=2.0d0*(q(2)*q(4)+q(1)*q(3)) 300 u(3,2)=2.0d0*(q(3)*q(4)-q(1)*q(2)) 301 u(3,3)=q(1)*q(1)-q(2)*q(2)-q(3)*q(3)+q(4)*q(4) 302c 303 rms1=zero 304 do 14 k=1,nx 305 do 15 i=1,3 306 q(i)=0.0d0 307 do 16 j=1,3 308 q(i)=q(i)+u(i,j)*(x(k,j)-xt(j)) 309 16 continue 310 15 continue 311 if(mod) then 312 do 17 i=1,3 313 x(k,i)=q(i)+yt(i) 314 17 continue 315 endif 316 sd=0.0d0 317 do 18 j=1,3 318 sd=sd+(q(j)-y(ix(k),j)+yt(j))**2 319 18 continue 320 sdev(ix(k))=sdev(ix(k))+sd 321 rms1=rms1+ws(ix(k))*sd 322 14 continue 323 if(lpara) call ga_dgop(mag_d11,rms1,1,'+') 324 rms1=sqrt(rms1/wnorms) 325c 326 if(mod.and.nw.gt.0) then 327 do 19 l=1,nw 328 do 20 k=1,na 329 do 21 i=1,3 330 q(i)=0.0d0 331 do 22 j=1,3 332 q(i)=q(i)+u(i,j)*(xw(l,k,j)-xt(j)) 333 22 continue 334 21 continue 335 do 23 i=1,3 336 xw(l,k,i)=q(i)+yt(i) 337 23 continue 338 20 continue 339 19 continue 340 endif 341c 342c write(*,'(a,6f12.6)') 'Super out',(x(1,i),i=1,3),(x(nx,i),i=1,3) 343 return 344 end 345 real*8 function angl(x,y,z) 346c 347 implicit none 348c 349 real*8 x(3),y(3),z(3) 350c 351 real*8 xy(3),zy(3),rxy,rzy,phi 352 integer i 353c 354 rxy=0.0d0 355 rzy=0.0d0 356 do 1 i=1,3 357 xy(i)=x(i)-y(i) 358 zy(i)=z(i)-y(i) 359 1 continue 360 rxy=xy(1)*xy(1)+xy(2)*xy(2)+xy(3)*xy(3) 361 rzy=zy(1)*zy(1)+zy(2)*zy(2)+zy(3)*zy(3) 362 phi=(xy(1)*zy(1)+xy(2)*zy(2)+xy(3)*zy(3))/sqrt(rxy*rzy) 363 if(phi.lt.-1.0d0) phi=-1.0d0 364 if(phi.gt.1.0d0) phi=1.0d0 365 angl=acos(phi) 366c 367 return 368 end 369 real*8 function atom_radius(number) 370c 371 implicit none 372 integer number 373c 374 real*8 radius(0:105) 375c 376 data radius / 99999.99, 377 + 0.35, 1.22, 1.23, 0.89, 0.88, 0.77, 0.70, 0.66, 0.58, 1.60, 378 + 1.40, 1.36, 1.25, 1.17, 1.10, 1.04, 0.99, 1.91, 2.03, 1.74, 379 + 1.44, 1.32, 1.22, 1.19, 1.17,1.165, 1.16, 1.15, 1.17, 1.25, 380 + 1.25, 1.22, 1.21, 1.17, 1.14, 1.98, 2.22, 1.92, 1.62, 1.45, 381 + 1.34, 1.29, 1.27, 1.24, 1.25, 1.28, 1.34, 1.41, 1.50, 1.40, 382 + 1.41, 1.37, 1.33, 2.09, 2.35, 1.98, 1.69, 1.65, 1.65, 1.64, 383 + 1.65, 1.66, 1.65, 1.61, 1.59, 1.59, 1.58, 1.57, 1.56, 1.56, 384 + 1.56, 1.44, 1.34, 1.30, 1.28, 1.26, 1.26, 1.29, 1.34, 1.44, 385 + 1.55, 1.54, 1.52, 1.53, 1.50, 2.20, 3.24, 2.68, 2.25, 2.16, 386 + 1.93, 1.66, 1.57, 1.81, 2.21, 1.43, 1.42, 1.40, 1.39, 1.38, 387 + 1.37, 1.36, 1.34, 1.30, 1.30 / 388c 389 atom_radius=1.1d-01*radius(number) 390c 391 return 392 end 393 subroutine povinc(lfn,xmin,xmax,ymin,ymax,zmin,zmax) 394c 395 implicit none 396c 397 integer lfn 398 real*8 xmin,xmax,ymin,ymax,zmin,zmax 399c 400 open(unit=lfn,file='camera.inc',form='formatted',status='new', 401 + err=9) 402 write(lfn,1000) 403 + 0.0d0,0.0d0,-1.0d1*max(abs(xmax),abs(ymax),abs(zmax), 404 + abs(xmin),abs(ymin),abs(zmin)), 405 + 5.0d-1*(xmax+xmin),5.0d-1*(ymax+ymin),5.0d-1*(zmax+zmin) 406c + 0.0d0, 0.0d0, 0.0d0, 407c + 1.3d0, 0.0d0, 0.0d0, 408c + 0.0d0, 0.0d0, 1.0d0, 409c + 0.0d0, 1.0d0, 0.0d0 410 1000 format('camera {',/, 411 + ' location <',f12.6,',',f12.6,',',f12.6,'>',/, 412 + ' look_at <',f12.6,',',f12.6,',',f12.6,'>',/, 413 + ' angle 20.0',/,'}') 414 write(lfn,1001) 0.0d1, 0.0d1,-5.0d1, 1.0d0,1.0d0,1.0d0 415 write(lfn,1001) -1.0d1, 2.0d1,-1.0d1, 1.0d0,1.0d0,1.0d0 416 write(lfn,1001) 1.0d1, 2.0d1,-1.0d1, 1.0d0,1.0d0,1.0d0 417 write(lfn,1001) 0.0d1, 1.0d1,-2.0d1, 1.0d0,1.0d0,1.0d0 418 1001 format('light_source { <',f12.6,',',f12.6,',',f12.6, 419 + '> color rgb <',f4.2,',',f4.2,',',f4.2,'> }') 420 write(lfn,1002) 0.0d0,0.0d0,0.0d0 421 1002 format('background { color rgb <',f4.2,',',f4.2,',',f4.2,'> }') 422 close(unit=lfn) 423c 424 9 continue 425 open(unit=lfn,file='colors.inc',form='formatted',status='new', 426 + err=99) 427 write(lfn,2000) 428 2000 format('#declare Colors_Inc_Temp = version ;',/,'#version 2.0 ;') 429 write(lfn,2001) 'Gray05',0.05,0.05,0.05 430 write(lfn,2001) 'Gray05',0.05,0.05,0.05 431 write(lfn,2001) 'Gray10',0.10,0.10,0.10 432 write(lfn,2001) 'Gray15',0.15,0.15,0.15 433 write(lfn,2001) 'Gray20',0.20,0.20,0.20 434 write(lfn,2001) 'Gray25',0.25,0.25,0.25 435 write(lfn,2001) 'Gray30',0.30,0.30,0.30 436 write(lfn,2001) 'Gray35',0.35,0.35,0.35 437 write(lfn,2001) 'Gray40',0.40,0.40,0.40 438 write(lfn,2001) 'Gray45',0.45,0.45,0.45 439 write(lfn,2001) 'Gray50',0.50,0.50,0.50 440 write(lfn,2001) 'Gray55',0.55,0.55,0.55 441 write(lfn,2001) 'Gray60',0.60,0.60,0.60 442 write(lfn,2001) 'Gray65',0.65,0.65,0.65 443 write(lfn,2001) 'Gray70',0.70,0.70,0.70 444 write(lfn,2001) 'Gray75',0.75,0.75,0.75 445 write(lfn,2001) 'Gray80',0.80,0.80,0.80 446 write(lfn,2001) 'Gray85',0.85,0.85,0.85 447 write(lfn,2001) 'Gray90',0.90,0.90,0.90 448 write(lfn,2001) 'Gray95',0.95,0.95,0.95 449 write(lfn,2001) 'DimGray',0.329412,0.329412,0.329412 450 write(lfn,2001) 'DimGrey',0.329412,0.329412,0.329412 451 write(lfn,2001) 'Gray',0.752941,0.752941,0.752941 452 write(lfn,2001) 'Grey',0.752941,0.752941,0.752941 453 write(lfn,2001) 'LightGray',0.658824,0.658824,0.658824 454 write(lfn,2001) 'LightGrey',0.658824,0.658824,0.658824 455 write(lfn,2001) 'VLightGrey',0.80,0.80,0.80 456 write(lfn,2001) 'White',1.0,1.0,1.0 457 write(lfn,2001) 'Red',1.0,0.0,0.0 458 write(lfn,2001) 'Green',0.0,1.0,0.0 459 write(lfn,2001) 'Blue',0.0,0.0,1.0 460 write(lfn,2001) 'Yellow',1.0,1.0,0.0 461 write(lfn,2001) 'Cyan',0.0,1.0,1.0 462 write(lfn,2001) 'Magenta',1.0,0.0,1.0 463 write(lfn,2001) 'Black',0.0,0.0,0.0 464 write(lfn,2001) 'Aquamarine',0.439216,0.858824,0.576471 465 write(lfn,2001) 'BlueViolet',0.62352,0.372549,0.623529 466 write(lfn,2001) 'Brown',0.647059,0.164706,0.164706 467 write(lfn,2001) 'CadetBlue',0.372549,0.623529,0.623529 468 write(lfn,2001) 'Coral',1.0,0.498039,0.0 469 write(lfn,2001) 'CornflowerBlue',0.258824,0.258824,0.435294 470 write(lfn,2001) 'DarkGreen',0.184314,0.309804,0.184314 471 write(lfn,2001) 'DarkOliveGreen',0.309804,0.309804,0.184314 472 write(lfn,2001) 'DarkOrchid',0.6,0.196078,0.8 473 write(lfn,2001) 'DarkSlateBlue',0.419608,0.137255,0.556863 474 write(lfn,2001) 'DarkSlateGray',0.184314,0.309804,0.309804 475 write(lfn,2001) 'DarkSlateGrey',0.184314,0.309804,0.309804 476 write(lfn,2001) 'DarkTurquoise',0.439216,0.576471,0.858824 477 write(lfn,2001) 'Firebrick',0.556863,0.137255,0.137255 478 write(lfn,2001) 'ForestGreen',0.137255,0.556863,0.137255 479 write(lfn,2001) 'Gold',0.8,0.498039,0.196078 480 write(lfn,2001) 'Goldenrod',0.858824,0.858824,0.439216 481 write(lfn,2001) 'GreenYellow',0.576471,0.858824,0.439216 482 write(lfn,2001) 'IndianRed',0.309804,0.184314,0.184314 483 write(lfn,2001) 'Khaki',0.623529,0.623529,0.372549 484 write(lfn,2001) 'LightBlue',0.74902,0.847059,0.847059 485 write(lfn,2001) 'LightSteelBlue',0.560784,0.560784,0.737255 486 write(lfn,2001) 'LimeGreen',0.196078,0.8,0.196078 487 write(lfn,2001) 'Maroon',0.556863,0.137255,0.419608 488 write(lfn,2001) 'MediumAquamarine',0.196078,0.8,0.6 489 write(lfn,2001) 'MediumBlue',0.196078,0.196078,0.8 490 write(lfn,2001) 'MediumForestGreen',0.419608,0.556863,0.137255 491 write(lfn,2001) 'MediumGoldenrod',0.917647,0.917647,0.678431 492 write(lfn,2001) 'MediumOrchid',0.576471,0.439216,0.858824 493 write(lfn,2001) 'MediumSeaGreen',0.258824,0.435294,0.258824 494 write(lfn,2001) 'MediumSlateBlue',0.498039,1.0,0.0 495 write(lfn,2001) 'MediumSpringGreen',0.498039,1.0,0.0 496 write(lfn,2001) 'MediumTurquoise',0.439216,0.858824,0.858824 497 write(lfn,2001) 'MediumVioletRed',0.858824,0.439216,0.576471 498 write(lfn,2001) 'MidnightBlue',0.184314,0.184314,0.309804 499 write(lfn,2001) 'Navy',0.137255,0.137255,0.556863 500 write(lfn,2001) 'NavyBlue',0.137255,0.137255,0.556863 501 write(lfn,2001) 'Orange',1,0.5,0.0 502 write(lfn,2001) 'OrangeRed',1.0,0.498039,0.0 503 write(lfn,2001) 'Orchid',0.858824,0.439216,0.858824 504 write(lfn,2001) 'PaleGreen',0.560784,0.737255,0.560784 505 write(lfn,2001) 'Pink',0.737255,0.560784,0.560784 506 write(lfn,2001) 'Plum',0.917647,0.678431,0.917647 507 write(lfn,2001) 'Salmon',0.435294,0.258824,0.258824 508 write(lfn,2001) 'SeaGreen',0.137255,0.556863,0.419608 509 write(lfn,2001) 'Sienna',0.556863,0.419608,0.137255 510 write(lfn,2001) 'SkyBlue',0.196078,0.6,0.8 511 write(lfn,2001) 'SlateBlue',0.0,0.498039,1.0 512 write(lfn,2001) 'SpringGreen',0.0,1.0,0.498039 513 write(lfn,2001) 'SteelBlue',0.137255,0.419608,0.556863 514 write(lfn,2001) 'Tan',0.858824,0.576471,0.439216 515 write(lfn,2001) 'Thistle',0.847059,0.74902,0.847059 516 write(lfn,2001) 'Turquoise',0.678431,0.917647,0.917647 517 write(lfn,2001) 'Violet',0.309804,0.184314,0.309804 518 write(lfn,2001) 'VioletRed',0.8,0.196078,0.6 519 write(lfn,2001) 'Wheat',0.847059,0.847059,0.74902 520 write(lfn,2001) 'YellowGreen',0.6,0.8,0.196078 521 write(lfn,2001) 'SummerSky',0.22,0.69,0.87 522 write(lfn,2001) 'RichBlue',0.35,0.35,0.67 523 write(lfn,2001) 'Brass',0.71,0.65,0.26 524 write(lfn,2001) 'Copper',0.72,0.45,0.20 525 write(lfn,2001) 'Bronze',0.55,0.47,0.14 526 write(lfn,2001) 'Bronze2',0.65,0.49,0.24 527 write(lfn,2001) 'Silver',0.90,0.91,0.98 528 write(lfn,2001) 'BrightGold',0.85,0.85,0.10 529 write(lfn,2001) 'OldGold',0.81,0.71,0.23 530 write(lfn,2001) 'Feldspar',0.82,0.57,0.46 531 write(lfn,2001) 'Quartz',0.85,0.85,0.95 532 write(lfn,2001) 'Mica',0.0,0.0,0.0 533 write(lfn,2001) 'NeonPink',1.00,0.43,0.78 534 write(lfn,2001) 'DarkPurple',0.53,0.12,0.47 535 write(lfn,2001) 'NeonBlue',0.30,0.30,1.00 536 write(lfn,2001) 'CoolCopper',0.85,0.53,0.10 537 write(lfn,2001) 'MandarinOrange',0.89,0.47,0.20 538 write(lfn,2001) 'LightWood',0.91,0.76,0.65 539 write(lfn,2001) 'MediumWood',0.65,0.50,0.39 540 write(lfn,2001) 'DarkWood',0.52,0.37,0.26 541 write(lfn,2001) 'SpicyPink',1.00,0.11,0.68 542 write(lfn,2001) 'SemiSweetChoc',0.42,0.26,0.15 543 write(lfn,2001) 'BakersChoc',0.36,0.20,0.09 544 write(lfn,2001) 'Flesh',0.96,0.80,0.69 545 write(lfn,2001) 'NewTan',0.92,0.78,0.62 546 write(lfn,2001) 'NewMidnightBlue',0.00,0.00,0.61 547 write(lfn,2001) 'VeryDarkBrown',0.35,0.16,0.14 548 write(lfn,2001) 'DarkBrown',0.36,0.25,0.20 549 write(lfn,2001) 'DarkTan',0.59,0.41,0.31 550 write(lfn,2001) 'GreenCopper',0.32,0.49,0.46 551 write(lfn,2001) 'DkGreenCopper',0.29,0.46,0.43 552 write(lfn,2001) 'DustyRose',0.52,0.39,0.39 553 write(lfn,2001) 'HuntersGreen',0.13,0.37,0.31 554 write(lfn,2001) 'Scarlet',0.55,0.09,0.09 555 write(lfn,2002) 'Clear',1.0,1.0,1.0,1.0 556 2001 format('#declare ',a,' = color red ',f8.6,' green ',f8.6, 557 + ' blue ',f8.6,' ;') 558 2002 format('#declare ',a,' = color red ',f8.6,' green ',f8.6, 559 + ' blue ',f8.6,' filter ',f8.6,' ;') 560 write(lfn,2003) 561 2003 format('#declare Plane_Map = 0 ;',/,'#declare Sphere_Map = 1 ;',/, 562 + '#declare Cylinder_Map = 2 ;',/,'#declare Torus_Map = 5 ;',/, 563 + '#declare Bi = 2 ;',/,'#declare Norm = 4 ;',/, 564 + '#version Colors_Inc_Temp ;') 565 close(unit=lfn) 566c 567 99 continue 568 open(unit=lfn,file='plane.inc',form='formatted',status='new', 569 + err=999) 570 write(lfn,3001) 571 3001 format('plane{ <0,1,0>,-1 pigment { White } }') 572 close(unit=lfn) 573c 574 999 continue 575c 576 return 577 end 578 real*8 function torsion(a,b,c,d) 579c 580 implicit none 581c 582 real*8 a(3),b(3),c(3),d(3) 583c 584 torsion=0.0d0 585c 586 return 587 end 588 character*255 function atom_color(number) 589 implicit none 590 integer number 591c 592 integer indexc(0:105) 593c 594 data indexc / 12, 595 + 7, 5, 14, 12, 13, 0, 1, 2, 6, 12, 596 + 13, 15, 9, 6, 8, 3, 13, 12, 16, 16, 597 + 12, 9, 12, 9, 9, 8, 12, 10, 10, 10, 598 + 12, 12, 12, 12, 10, 12, 12, 12, 12, 12, 599 + 12, 12, 12, 12, 12, 12, 16, 12, 12, 12, 600 + 12, 12, 11, 12, 12, 8, 12, 12, 12, 12, 601 + 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 602 + 12, 12, 12, 12, 12, 12, 12, 16, 17, 16, 603 + 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 604 + 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 605 + 12, 12, 12, 12, 12 / 606c 607 atom_color='white ' 608 if(indexc(number).eq.0) atom_color='LightGrey ' 609 if(indexc(number).eq.1) atom_color='SkyBlue ' 610 if(indexc(number).eq.2) atom_color='Red ' 611 if(indexc(number).eq.3) atom_color='Yellow ' 612 if(indexc(number).eq.4) atom_color='White ' 613 if(indexc(number).eq.5) atom_color='Pink ' 614 if(indexc(number).eq.6) atom_color='Goldenrod ' 615 if(indexc(number).eq.7) atom_color='Blue ' 616 if(indexc(number).eq.8) atom_color='Orange ' 617 if(indexc(number).eq.9) atom_color='Gray85 ' 618 if(indexc(number).eq.10) atom_color='Brown ' 619 if(indexc(number).eq.11) atom_color='Purple ' 620 if(indexc(number).eq.12) atom_color='SpicyPink ' 621 if(indexc(number).eq.13) atom_color='Green ' 622 if(indexc(number).eq.14) atom_color='Firebrick ' 623 if(indexc(number).eq.15) atom_color='DarkGreen ' 624 if(indexc(number).eq.16) atom_color='Silver ' 625 if(indexc(number).eq.17) atom_color='Gold ' 626c 627 return 628 end 629 subroutine md_jacobi(a,n,na,d,v,nrot) 630c 631c compute eigenvectors and eigenvalues for real symmetric 632c matrix using the Jacobi diagonalization 633c 634 implicit none 635c 636 integer nrmax 637 parameter(nrmax=100) 638c 639 real*8 zero,half,one,two 640 parameter(zero=0.0d0) 641 parameter(half=0.5d0) 642 parameter(one=1.0d0) 643 parameter(two=2.0d0) 644c 645 integer n,na,nrot 646 real*8 a(na,na),d(na),v(na,na) 647 real*8 at,b,dma,q 648c 649 integer i,j,k,l 650 real*8 c,s,t,sum,temp 651c 652 do 1 i=1,n 653 do 2 j=1,n 654 v(i,j)=zero 655 2 continue 656 v(i,i)=one 657 d(i)=a(i,i) 658 1 continue 659c 660 nrot=0 661 do 3 l=1,nrmax 662 nrot=nrot+1 663 sum=zero 664 do 4 i=1,n-1 665 do 5 j=i+1,n 666 sum=sum+abs(a(i,j)) 667 5 continue 668 4 continue 669 if(sum.eq.zero) then 670 do 6 i=1,n-1 671 do 7 j=i+1,n 672 if(d(i).gt.d(j)) then 673 temp=d(i) 674 d(i)=d(j) 675 d(j)=temp 676 do 8 k=1,n 677 temp=v(k,i) 678 v(k,i)=v(k,j) 679 v(k,j)=temp 680 8 continue 681 endif 682 7 continue 683 6 continue 684 return 685 endif 686 do 9 j=2,n 687 do 10 i=1,j-1 688 b=a(i,j) 689 if(abs(b).gt.zero) then 690 dma=d(j)-d(i) 691 if(abs(dma)+abs(b).le.abs(dma)) then 692 t=b/dma 693 else 694 q=half*dma/b 695 t=sign(one/(abs(q)+sqrt(one+q*q)),q) 696 endif 697 c=one/sqrt(t*t+one) 698 s=t*c 699 a(i,j)=zero 700 do 11 k=1,i-1 701 at=c*a(k,i)-s*a(k,j) 702 a(k,j)=s*a(k,i)+c*a(k,j) 703 a(k,i)=at 704 11 continue 705 do 12 k=i+1,j-1 706 at=c*a(i,k)-s*a(k,j) 707 a(k,j)=s*a(i,k)+c*a(k,j) 708 a(i,k)=at 709 12 continue 710 do 13 k=j+1,n 711 at=c*a(i,k)-s*a(j,k) 712 a(j,k)=s*a(i,k)+c*a(j,k) 713 a(i,k)=at 714 13 continue 715 do 14 k=1,n 716 at=c*v(k,i)-s*v(k,j) 717 v(k,j)=s*v(k,i)+c*v(k,j) 718 v(k,i)=at 719 14 continue 720 at=c*c*d(i)+s*s*d(j)-two*c*s*b 721 d(j)=s*s*d(i)+c*c*d(j)+two*c*s*b 722 d(i)=at 723 endif 724 10 continue 725 9 continue 726 3 continue 727c 728 call md_abort('md_jacobi: maximum iterations reached',0) 729c 730 return 731 end 732 subroutine swatch(today,now) 733c 734 implicit none 735c 736 character*10 today,now 737c 738#if defined(LINUX) 739 character*26 string 740#endif 741#if defined(IBM) 742 character*26 string 743#endif 744#if defined(SP1) || defined(SOLARIS) 745 character*26 string 746#endif 747#if defined(SGI) 748 character*9 string 749#endif 750c 751 today='00/00/00 ' 752 now='00:00:00' 753c 754#if defined(IBM) 755 call fdate(string) 756 if(string(4:6).eq.'Jan') today(1:2)='01' 757 if(string(4:6).eq.'Feb') today(1:2)='02' 758 if(string(4:6).eq.'Mar') today(1:2)='03' 759 if(string(4:6).eq.'Apr') today(1:2)='04' 760 if(string(4:6).eq.'May') today(1:2)='05' 761 if(string(4:6).eq.'Jun') today(1:2)='06' 762 if(string(4:6).eq.'Jul') today(1:2)='07' 763 if(string(4:6).eq.'Aug') today(1:2)='08' 764 if(string(4:6).eq.'Sep') today(1:2)='09' 765 if(string(4:6).eq.'Oct') today(1:2)='10' 766 if(string(4:6).eq.'Nov') today(1:2)='11' 767 if(string(4:6).eq.'Dec') today(1:2)='12' 768 today(7:8)=string(8:9) 769 today(4:5)=string(1:2) 770 now=string(11:20) 771#endif 772#if defined(SP1) || defined(SOLARIS) 773 call util_date(string) 774 if(string(5:7).eq.'Jan') today(1:2)='01' 775 if(string(5:7).eq.'Feb') today(1:2)='02' 776 if(string(5:7).eq.'Mar') today(1:2)='03' 777 if(string(5:7).eq.'Apr') today(1:2)='04' 778 if(string(5:7).eq.'May') today(1:2)='05' 779 if(string(5:7).eq.'Jun') today(1:2)='06' 780 if(string(5:7).eq.'Jul') today(1:2)='07' 781 if(string(5:7).eq.'Aug') today(1:2)='08' 782 if(string(5:7).eq.'Sep') today(1:2)='09' 783 if(string(5:7).eq.'Oct') today(1:2)='10' 784 if(string(5:7).eq.'Nov') today(1:2)='11' 785 if(string(5:7).eq.'Dec') today(1:2)='12' 786 today(7:8)=string(23:24) 787 today(4:5)=string(9:10) 788 now=string(11:20) 789#endif 790#if defined(LINUX) 791 call util_date(string) 792 if(string(5:7).eq.'Jan') today(1:2)='01' 793 if(string(5:7).eq.'Feb') today(1:2)='02' 794 if(string(5:7).eq.'Mar') today(1:2)='03' 795 if(string(5:7).eq.'Apr') today(1:2)='04' 796 if(string(5:7).eq.'May') today(1:2)='05' 797 if(string(5:7).eq.'Jun') today(1:2)='06' 798 if(string(5:7).eq.'Jul') today(1:2)='07' 799 if(string(5:7).eq.'Aug') today(1:2)='08' 800 if(string(5:7).eq.'Sep') today(1:2)='09' 801 if(string(5:7).eq.'Oct') today(1:2)='10' 802 if(string(5:7).eq.'Nov') today(1:2)='11' 803 if(string(5:7).eq.'Dec') today(1:2)='12' 804 today(7:8)=string(23:24) 805 today(4:5)=string(9:10) 806 now=string(11:20) 807#endif 808#if defined(SGI) 809 call date(string) 810 if(string(4:6).eq.'Jan') today(1:2)='01' 811 if(string(4:6).eq.'Feb') today(1:2)='02' 812 if(string(4:6).eq.'Mar') today(1:2)='03' 813 if(string(4:6).eq.'Apr') today(1:2)='04' 814 if(string(4:6).eq.'May') today(1:2)='05' 815 if(string(4:6).eq.'Jun') today(1:2)='06' 816 if(string(4:6).eq.'Jul') today(1:2)='07' 817 if(string(4:6).eq.'Aug') today(1:2)='08' 818 if(string(4:6).eq.'Sep') today(1:2)='09' 819 if(string(4:6).eq.'Oct') today(1:2)='10' 820 if(string(4:6).eq.'Nov') today(1:2)='11' 821 if(string(4:6).eq.'Dec') today(1:2)='12' 822 today(7:8)=string(8:9) 823 today(4:5)=string(1:2) 824 call time(now(1:8)) 825 now(9:10)=' ' 826#endif 827 if(today(4:4).eq.' ') today(4:4)='0' 828 return 829 end 830 subroutine matinv(a,n,ndim) 831c 832 implicit none 833 834 integer maxdim 835 real*8 zero,small,one 836 parameter(maxdim=3) 837 parameter(zero=0.0d0) 838 parameter(small=1.0d-6) 839 parameter(one=1.0d0) 840c 841 integer n,ndim 842 real*8 a(ndim,ndim) 843 integer ia(2,maxdim),ib(maxdim),ic(maxdim) 844 real*8 d(maxdim) 845c 846 integer idim,i,j,k,l,m 847 real*8 b,e 848c 849 if(ndim.gt.maxdim) call md_abort('matinv dimension error',0) 850c 851 do 1 idim=1,n 852 ia(1,idim)=0 853 ia(2,idim)=0 854 1 continue 855c 856 do 9 idim=1,n 857 b=zero 858 do 3 l=1,n 859 do 4 m=1,n 860 if(ia(1,l).ne.1.and.ia(2,m).ne.1) then 861 e=dabs(a(l,m)) 862 if(e.ge.b) then 863 i=l 864 k=m 865 endif 866 8 b=dmax1(b,e) 867 endif 868 4 continue 869 3 continue 870 ia(1,i)=1 871 ia(2,k)=1 872 ib(k)=i 873 ic(i)=k 874 b=a(i,k) 875c 876 if(dabs(b).lt.small) call md_abort('arg_matinv singular matrix',0) 877 a(i,k)=one/b 878 do 6 l=1,n 879 if(l.ne.k) a(i,l)=-a(i,l)/b 880 6 continue 881 do 5 l=1,n 882 do 7 m=1,n 883 if(l.ne.i.and.m.ne.k) a(l,m)=a(l,m)+a(l,k)*a(i,m) 884 7 continue 885 5 continue 886 do 11 l=1,n 887 if(l.ne.i) a(l,k)=a(l,k)/b 888 11 continue 889 9 continue 890c 891 do 15 l=1,n 892 do 13 j=1,n 893 k=ib(j) 894 d(j)=a(k,l) 895 13 continue 896 do 14 j=1,n 897 a(j,l)=d(j) 898 14 continue 899 15 continue 900c 901 do 16 l=1,n 902 do 17 j=1,n 903 k=ic(j) 904 d(j)=a(l,k) 905 17 continue 906 do 18 j=1,n 907 a(l,j)=d(j) 908 18 continue 909 16 continue 910c 911 return 912 end 913 logical function frequency(istep,nstep) 914c 915 implicit none 916c 917 integer istep,nstep 918c 919 if(nstep.le.0) then 920 frequency=.false. 921 else 922 frequency=mod(istep,nstep).eq.0 923 endif 924c 925 return 926 end 927 subroutine rolex(elaps,cputim) 928c 929 implicit none 930c 931 real*8 elaps,cputim 932c 933 real*8 util_wallsec,util_cpusec 934 external util_wallsec,util_cpusec 935c 936 elaps=util_wallsec() 937 cputim=util_cpusec() 938c 939 return 940 end 941 subroutine timer_init() 942c 943 implicit none 944c 945 call timer(0,0) 946c 947 return 948 end 949 subroutine timer_reset(itime) 950c 951 implicit none 952c 953 integer itime 954c 955 call timer(itime,-2) 956c 957 return 958 end 959 subroutine timer_start(itime) 960c 961 implicit none 962c 963 integer itime 964c 965 call timer(itime,0) 966c 967 return 968 end 969 subroutine timer_stop(itime) 970c 971 implicit none 972c 973 integer itime 974c 975 call timer(itime,1) 976c 977 return 978 end 979 subroutine timer(itime,iopt) 980c 981 implicit none 982c 983 integer itime,iopt 984c 985 real*8 elaps,cputim 986c 987 integer mtime 988 parameter(mtime=250) 989 integer ncall(250) 990 real*8 ttime(250,3),ctime(250,3) 991 common/tim/ncall,ttime,ctime 992c 993 integer i 994c 995 call rolex(elaps,cputim) 996c 997 if(itime.eq.0) then 998 do 1 i=1,mtime 999 ncall(i)=0 1000 ctime(i,1)=0.0d0 1001 ttime(i,1)=0.0d0 1002 ctime(i,2)=0.0d0 1003 ttime(i,2)=0.0d0 1004 ctime(i,3)=1.0d9 1005 ttime(i,3)=1.0d9 1006 1 continue 1007 elseif(itime.le.0.or.itime.gt.mtime) then 1008 call md_abort('Timer index out of range',0) 1009 elseif(iopt.eq.-2) then 1010 ncall(itime)=0 1011 ctime(itime,1)=0.0d0 1012 ttime(itime,1)=0.0d0 1013 ctime(itime,2)=0.0d0 1014 ttime(itime,2)=0.0d0 1015 ctime(itime,3)=1.0d9 1016 ttime(itime,3)=1.0d9 1017 elseif(iopt.eq.-1) then 1018 ncall(itime)=0 1019 ctime(itime,1)=-cputim 1020 ttime(itime,1)=-elaps 1021 ctime(itime,2)=-cputim 1022 ttime(itime,2)=-elaps 1023 ctime(itime,3)=1.0d9 1024 ttime(itime,3)=1.0d9 1025 elseif(iopt.eq.0) then 1026 ctime(itime,1)=ctime(itime,1)-cputim 1027 ttime(itime,1)=ttime(itime,1)-elaps 1028 ctime(itime,2)=-cputim 1029 ttime(itime,2)=-elaps 1030 elseif(iopt.eq.1) then 1031 ncall(itime)=ncall(itime)+1 1032 ctime(itime,1)=ctime(itime,1)+cputim 1033 ttime(itime,1)=ttime(itime,1)+elaps 1034 ctime(itime,2)=ctime(itime,2)+cputim 1035 ttime(itime,2)=ttime(itime,2)+elaps 1036 ctime(itime,3)=min(ctime(itime,2),ctime(itime,3)) 1037 ttime(itime,3)=min(ttime(itime,2),ttime(itime,3)) 1038 else 1039 call md_abort('Unimplemented timer option',0) 1040 endif 1041c 1042 return 1043 end 1044 real*8 function timer_cpu(itime) 1045c 1046 implicit none 1047c 1048 integer itime 1049c 1050 integer mtime 1051 parameter(mtime=250) 1052 integer ncall(250) 1053 real*8 ttime(250,3),ctime(250,3) 1054 common/tim/ncall,ttime,ctime 1055c 1056 if(itime.le.0.or.itime.gt.mtime) 1057 + call md_abort('Illegal timer index',0) 1058c 1059 if(ncall(itime).le.0) then 1060 timer_cpu=0.0d0 1061 else 1062 timer_cpu=ctime(itime,2) 1063 endif 1064c 1065 return 1066 end 1067 real*8 function timer_wall(itime) 1068c 1069 implicit none 1070c 1071 integer itime 1072c 1073 integer mtime 1074 parameter(mtime=250) 1075 integer ncall(250) 1076 real*8 ttime(250,3),ctime(250,3) 1077 common/tim/ncall,ttime,ctime 1078c 1079 if(itime.le.0.or.itime.gt.mtime) 1080 + call md_abort('Illegal timer index',0) 1081c 1082 if(ncall(itime).le.0) then 1083 timer_wall=0.0d0 1084 else 1085 timer_wall=ttime(itime,2) 1086 endif 1087c 1088 return 1089 end 1090 integer function timer_calls(itime) 1091c 1092 integer itime 1093c 1094 integer mtime 1095 parameter(mtime=250) 1096 integer ncall(250) 1097 real*8 ttime(250,3),ctime(250,3) 1098 common/tim/ncall,ttime,ctime 1099c 1100 if(itime.le.0.or.itime.gt.mtime) 1101 + call md_abort('Illegal timer index',0) 1102c 1103 timer_calls=ncall(itime) 1104c 1105 return 1106 end 1107 real*8 function timer_cpu_minimum(itime) 1108c 1109 implicit none 1110c 1111 integer itime 1112c 1113 integer mtime 1114 parameter(mtime=250) 1115 integer ncall(250) 1116 real*8 ttime(250,3),ctime(250,3) 1117 common/tim/ncall,ttime,ctime 1118c 1119 if(itime.le.0.or.itime.gt.mtime) 1120 + call md_abort('Illegal timer index',0) 1121c 1122 if(ncall(itime).le.0) then 1123 timer_cpu_minimum=0.0d0 1124 else 1125 timer_cpu_minimum=ctime(itime,3) 1126 endif 1127c 1128 return 1129 end 1130 real*8 function timer_wall_minimum(itime) 1131c 1132 implicit none 1133c 1134 integer itime 1135c 1136 integer mtime 1137 parameter(mtime=250) 1138 integer ncall(250) 1139 real*8 ttime(250,3),ctime(250,3) 1140 common/tim/ncall,ttime,ctime 1141c 1142 if(itime.le.0.or.itime.gt.mtime) 1143 + call md_abort('Illegal timer index',0) 1144c 1145 if(ncall(itime).le.0) then 1146 timer_wall_minimum=0.0d0 1147 else 1148 timer_wall_minimum=ttime(itime,3) 1149 endif 1150c 1151 return 1152 end 1153 real*8 function timer_cpu_average(itime) 1154c 1155 implicit none 1156c 1157 integer itime 1158c 1159 integer mtime 1160 parameter(mtime=250) 1161 integer ncall(250) 1162 real*8 ttime(250,3),ctime(250,3) 1163 common/tim/ncall,ttime,ctime 1164c 1165 if(itime.le.0.or.itime.gt.mtime) 1166 + call md_abort('Illegal timer index',0) 1167c 1168 if(ncall(itime).le.0) then 1169 timer_cpu_average=0.0d0 1170 else 1171 timer_cpu_average=ctime(itime,1)/dble(ncall(itime)) 1172 endif 1173c 1174 return 1175 end 1176 real*8 function timer_wall_average(itime) 1177c 1178 implicit none 1179c 1180 integer itime 1181c 1182 integer mtime 1183 parameter(mtime=250) 1184 integer ncall(250) 1185 real*8 ttime(250,3),ctime(250,3) 1186 common/tim/ncall,ttime,ctime 1187c 1188 if(itime.le.0.or.itime.gt.mtime) 1189 + call md_abort('Illegal timer index',0) 1190c 1191 if(ncall(itime).le.0) then 1192 timer_wall_average=0.0d0 1193 else 1194 timer_wall_average=ttime(itime,1)/dble(ncall(itime)) 1195 endif 1196c 1197 return 1198 end 1199 real*8 function timer_cpu_total(itime) 1200c 1201 implicit none 1202c 1203 integer itime 1204c 1205 integer mtime 1206 parameter(mtime=250) 1207 integer ncall(250) 1208 real*8 ttime(250,3),ctime(250,3) 1209 common/tim/ncall,ttime,ctime 1210c 1211 if(itime.le.0.or.itime.gt.mtime) 1212 + call md_abort('Illegal timer index',0) 1213c 1214 if(ncall(itime).le.0) then 1215 timer_cpu_total=0.0d0 1216 else 1217 timer_cpu_total=ctime(itime,1) 1218 endif 1219c 1220 return 1221 end 1222 real*8 function timer_wall_total(itime) 1223c 1224 implicit none 1225c 1226 integer itime 1227c 1228 integer mtime 1229 parameter(mtime=250) 1230 integer ncall(250) 1231 real*8 ttime(250,3),ctime(250,3) 1232 common/tim/ncall,ttime,ctime 1233c 1234 if(itime.le.0.or.itime.gt.mtime) 1235 + call md_abort('Illegal timer index',0) 1236c 1237 if(ncall(itime).le.0) then 1238 timer_wall_total=0.0d0 1239 else 1240 timer_wall_total=ttime(itime,1) 1241 endif 1242c 1243 return 1244 end 1245 subroutine error(lauto,lapprox,maxacf,data,ndata, 1246 + aver,drift,stderr,corerr,ratio) 1247c 1248 implicit none 1249c 1250 real*8 zero,one 1251 parameter(zero=0.0d0) 1252 parameter(one=1.0d0) 1253c 1254#include "mafdecls.fh" 1255c 1256 logical lauto,lapprox 1257 integer ndata,lenacf,maxacf 1258 real*8 data(ndata),aver,drift,stderr,corerr,ratio 1259c 1260 integer i,i_acf,l_acf 1261 real*8 dsum,ddsum,dtsum,tsum,ttsum,dstep 1262c 1263c integer nacf,kapprx(15),iapprx,klarge 1264c real*8 warg 1265c real*8 data(ndata),acf(ndata),approx(15),cdac(15),weight 1266c 1267c integer i,j,k,l,m,nacfa,nfunc 1268c real*8 dsum,ddsum,dtsum,tsum,ttsum,dstep,dfsum,dvar 1269c real*8 cdap(15,15),cdaq(15,15),cdad(15),rfact 1270c real*8 xappm,xappi,xappj,wterm,wsum1,wsum2,cdawgt 1271c 1272c arg_error iopt = 0 : average 1273c drift 1274c standard error 1275c = 1 : average 1276c drift 1277c standard error 1278c autocorrelation function 1279c correlation error from actual acf 1280c sampling ration from actual acf 1281c = 2 : average 1282c drift 1283c standard error 1284c autocorrelation function 1285c correlation error from approximated acf 1286c sampling ration from approximated acf 1287c data : 1-dimensional array with data 1288c ndata : number of data 1289c nacf : length of autocorrelation function 1290c 1291c aver : average 1292c stderr : standard error 1293c drift : drift 1294c acf : autocorrelation function 1295c corerr : correlated error 1296c ratio : sampling ratio 1297c 1298 lenacf=maxacf 1299c 1300 dsum=zero 1301 ddsum=zero 1302 dtsum=zero 1303 tsum=zero 1304 ttsum=zero 1305 do 1 i=1,ndata 1306 dstep=dble(i) 1307 dsum=dsum+data(i) 1308 ddsum=ddsum+data(i)*data(i) 1309 dtsum=dtsum+dstep*data(i) 1310 tsum=tsum+dstep 1311 ttsum=ttsum+dstep*dstep 1312 1 continue 1313c 1314c average, drift and standard error 1315c 1316 aver=dsum/dble(ndata) 1317 drift=(dble(ndata)*dtsum-dsum*tsum)/(dble(ndata)*ttsum-tsum*tsum) 1318 stderr=sqrt(abs((ddsum/dble(ndata)-dsum*dsum/dble(ndata*ndata)))/ 1319 + (dble(ndata-1))) 1320c 1321 corerr=stderr 1322 ratio=one 1323c 1324 if(.not.lauto) return 1325c 1326 if(.not.ma_push_get(mt_dbl,lenacf,'acf',l_acf,i_acf)) 1327 + call md_abort('Failed to allocate memory for acf',0) 1328c 1329 call auto_corr(data,ndata,aver,dbl_mb(i_acf),lenacf,ratio) 1330c 1331 if(.not.ma_pop_stack(l_acf)) 1332 + call md_abort('Failed to deallocate memory for acf',0) 1333c 1334 corerr=ratio*stderr 1335c 1336 return 1337 end 1338 subroutine auto_corr(data,ndata,aver,acf,lacf,ratio) 1339c 1340 implicit none 1341c 1342 real*8 zero,half,one,two 1343 parameter(zero=0.0d0) 1344 parameter(half=5.0d-1) 1345 parameter(one=1.0d0) 1346 parameter(two=2.0d0) 1347c 1348 integer ndata,lacf,nacf 1349 real*8 data(ndata),acf(lacf),aver,ratio 1350c 1351 integer i,j 1352 real*8 dsum,dvar 1353c 1354 nacf=min(ndata,lacf) 1355c 1356 dsum=zero 1357 do 1 i=1,ndata 1358 dsum=dsum+(data(i)-aver)**2 1359 1 continue 1360 dvar=dble(ndata)/dsum 1361 do 2 i=1,nacf-1 1362 dsum=zero 1363 do 3 j=1,ndata-i 1364 dsum=dsum+(data(j)-aver)*(data(i+j)-aver) 1365 3 continue 1366 acf(i)=dvar*(dsum/dble(ndata-i))*(one-dble(i-1)/dble(ndata-2)) 1367 2 continue 1368c 1369 ratio=half 1370 do 4 i=1,nacf-1 1371 ratio=ratio+(one-(dble(i)/dble(ndata))**2)*abs(acf(i)) 1372 4 continue 1373 ratio=sqrt(two*abs(ratio)) 1374 if(ratio.lt.one) ratio=one 1375c 1376 lacf=nacf 1377c 1378 return 1379 end 1380 subroutine acf_approx(acf,acfapp,lacf) 1381c 1382 implicit none 1383c 1384 integer lacf 1385 real*8 acf(lacf),acfapp(lacf) 1386c 1387c integer kapprox(15) 1388c 1389c rfact=two*sqrt(dble(klarge))/dble(nacfa-1) 1390c wsum1=zero 1391c wsum2=zero 1392c warg=zero 1393c do 5 i=1,nacfa-1 1394c acf(i)=abs(acf(i)) 1395c if(acf(i).gt.zero) then 1396c wterm=exp(weight*dble(nacfa-i)/dble(nacfa-1)) 1397c wsum1=wsum1+wterm 1398c wsum2=wsum2+wterm*log(acf(i))/(rfact*dble(i)) 1399c endif 1400c 5 continue 1401c if(abs(wsum1).gt.small) warg=wsum2/wsum1 1402c do 55 i=1,nacfa-1 1403c acf(i)=acf(i)-exp(warg*dble(i)*rfact) 1404c 55 continue 1405c nfunc=iapprx 1406c if(nfunc.le.0) nfunc=15 1407c do 6 k=1,nfunc 1408c do 7 l=1,nfunc 1409c cdap(k,l)=zero 1410c do 8 m=1,nacfa-1 1411c xappm=dble(m)*rfact 1412c cdawgt=exp(weight*dble(nacfa-m)/dble(nacfa-1)) 1413c cdap(k,l)=cdap(k,l)+cdawgt*exp((-xappm)*xappm)* 1414c + approx(k)*approx(l)*(xappm**(kapprx(k)+kapprx(l))) 1415c 8 continue 1416c cdaq(k,l)=cdap(k,l) 1417c 7 continue 1418c 6 continue 1419c do 9 i=1,nfunc 1420c cdad(i)=zero 1421c do 10 j=1,nacfa-1 1422c xappj=dble(j)*rfact 1423c cdawgt=exp(weight*dble(nacfa-j)/dble(nacfa-1)) 1424c cdad(i)=cdad(i)+cdawgt*exp((-half)*xappj*xappj)* 1425c + acf(j)*approx(i)*(xappj**kapprx(i)) 1426c 10 continue 1427c 9 continue 1428c do 11 k=1,nfunc 1429c do 12 i=k,nfunc 1430c cdaq(i,k)=cdap(i,k) 1431c do 13 l=1,k-1 1432c cdaq(i,k)=cdaq(i,k)-cdaq(i,l)*cdaq(l,k) 1433c 13 continue 1434c 12 continue 1435c do 14 i=k+1,nfunc 1436c cdaq(k,i)=cdap(k,i) 1437c do 15 l=1,k-1 1438c cdaq(k,i)=cdaq(k,i)-cdaq(k,l)*cdaq(l,i) 1439c 15 continue 1440c cdaq(k,i)=cdaq(k,i)/cdaq(k,k) 1441c 14 continue 1442c 11 continue 1443c do 16 j=1,nfunc 1444c cdac(j)=cdad(j) 1445c do 17 i=1,j-1 1446c cdac(j)=cdac(j)-cdaq(j,i)*cdac(i) 1447c 17 continue 1448c cdac(j)=cdac(j)/cdaq(j,j) 1449c 16 continue 1450c do 18 k=1,nfunc 1451c j=nfunc+1-k 1452c do 19 i=j+1,nfunc 1453c cdac(j)=cdac(j)-cdaq(j,i)*cdac(i) 1454c 19 continue 1455c 18 continue 1456c do 20 i=1,nacfa-1 1457c xappi=dble(i)*rfact 1458c acf(i)=exp(warg*xappi) 1459c do 21 j=1,nfunc 1460c acf(i)=acf(i)+exp((-half)*xappi*xappi)* 1461c + approx(j)*cdac(j)*(xappi**kapprx(j)) 1462c 21 continue 1463c 20 continue 1464c 1465c autocorrelation function upto lag min(nacf,ndata) 1466c 1467c if(iopt.gt.0) then 1468c warg=zero 1469c nacfa=nacf 1470c if(nacfa.gt.ndata) nacfa=ndata 1471c dfsum=zero 1472c do 2 i=1,ndata 1473c dfsum=dfsum+(data(i)-aver)**2 1474c 2 continue 1475c dvar=dble(ndata)/dfsum 1476c do 3 i=1,nacfa-1 1477c dfsum=zero 1478c do 4 j=1,ndata-i 1479c dfsum=dfsum+(data(j)-aver)*(data(i+j)-aver) 1480c 4 continue 1481c acf(i)=dvar*(dfsum/dble(ndata-i))*(one-dble(i-1)/dble(ndata-2)) 1482c 3 continue 1483c endif 1484cc 1485cc approximate acf here 1486cc 1487c if(iopt.gt.1) then 1488c rfact=two*sqrt(dble(klarge))/dble(nacfa-1) 1489c wsum1=zero 1490c wsum2=zero 1491c warg=zero 1492c do 5 i=1,nacfa-1 1493c acf(i)=abs(acf(i)) 1494c if(acf(i).gt.zero) then 1495c wterm=exp(weight*dble(nacfa-i)/dble(nacfa-1)) 1496c wsum1=wsum1+wterm 1497c wsum2=wsum2+wterm*log(acf(i))/(rfact*dble(i)) 1498c endif 1499c 5 continue 1500c if(abs(wsum1).gt.small) warg=wsum2/wsum1 1501c do 55 i=1,nacfa-1 1502c acf(i)=acf(i)-exp(warg*dble(i)*rfact) 1503c 55 continue 1504c nfunc=iapprx 1505c if(nfunc.le.0) nfunc=15 1506c do 6 k=1,nfunc 1507c do 7 l=1,nfunc 1508c cdap(k,l)=zero 1509c do 8 m=1,nacfa-1 1510c xappm=dble(m)*rfact 1511c cdawgt=exp(weight*dble(nacfa-m)/dble(nacfa-1)) 1512c cdap(k,l)=cdap(k,l)+cdawgt*exp((-xappm)*xappm)* 1513c + approx(k)*approx(l)*(xappm**(kapprx(k)+kapprx(l))) 1514c 8 continue 1515c cdaq(k,l)=cdap(k,l) 1516c 7 continue 1517c 6 continue 1518c do 9 i=1,nfunc 1519c cdad(i)=zero 1520c do 10 j=1,nacfa-1 1521c xappj=dble(j)*rfact 1522c cdawgt=exp(weight*dble(nacfa-j)/dble(nacfa-1)) 1523c cdad(i)=cdad(i)+cdawgt*exp((-half)*xappj*xappj)* 1524c + acf(j)*approx(i)*(xappj**kapprx(i)) 1525c 10 continue 1526c 9 continue 1527c do 11 k=1,nfunc 1528c do 12 i=k,nfunc 1529c cdaq(i,k)=cdap(i,k) 1530c do 13 l=1,k-1 1531c cdaq(i,k)=cdaq(i,k)-cdaq(i,l)*cdaq(l,k) 1532c 13 continue 1533c 12 continue 1534c do 14 i=k+1,nfunc 1535c cdaq(k,i)=cdap(k,i) 1536c do 15 l=1,k-1 1537c cdaq(k,i)=cdaq(k,i)-cdaq(k,l)*cdaq(l,i) 1538c 15 continue 1539c cdaq(k,i)=cdaq(k,i)/cdaq(k,k) 1540c 14 continue 1541c 11 continue 1542c do 16 j=1,nfunc 1543c cdac(j)=cdad(j) 1544c do 17 i=1,j-1 1545c cdac(j)=cdac(j)-cdaq(j,i)*cdac(i) 1546c 17 continue 1547c cdac(j)=cdac(j)/cdaq(j,j) 1548c 16 continue 1549c do 18 k=1,nfunc 1550c j=nfunc+1-k 1551c do 19 i=j+1,nfunc 1552c cdac(j)=cdac(j)-cdaq(j,i)*cdac(i) 1553c 19 continue 1554c 18 continue 1555c do 20 i=1,nacfa-1 1556c xappi=dble(i)*rfact 1557c acf(i)=exp(warg*xappi) 1558c do 21 j=1,nfunc 1559c acf(i)=acf(i)+exp((-half)*xappi*xappi)* 1560c + approx(j)*cdac(j)*(xappi**kapprx(j)) 1561c 21 continue 1562c 20 continue 1563c endif 1564cc 1565cc sampling ratio 1566cc 1567c ratio=one 1568c if(iopt.gt.0) then 1569c ratio=half 1570c do 22 i=1,nacfa-1 1571c ratio=ratio+(one-(dble(i)/dble(ndata))**2)*abs(acf(i)) 1572c 22 continue 1573c ratio=sqrt(two*abs(ratio)) 1574c if(ratio.lt.one) ratio=one 1575c endif 1576cc 1577cc 1578c corerr=ratio*stderr 1579cc 1580 return 1581 end 1582 logical function md_zmat(x1,x2,x3,x4,d,a,t) 1583c 1584 implicit none 1585c 1586 real*8 x1(3),x2(3),x3(3),x4(3),d,a,t 1587c 1588 real*8 p(3),q(3),s(3),v(3),w(3) 1589 real*8 r 1590 integer i 1591c 1592 do 1 i=1,3 1593 p(i)=x3(i)-x2(i) 1594 q(i)=x4(i)-x2(i) 1595 1 continue 1596 v(1)=p(2)*q(3)-p(3)*q(2)+x2(1) 1597 v(2)=p(3)*q(1)-p(1)*q(3)+x2(2) 1598 v(3)=p(1)*q(2)-p(2)*q(1)+x2(3) 1599c 1600 r=sqrt(p(1)*p(1)+p(2)*p(2)+p(3)*p(3)) 1601 do 2 i=1,3 1602 s(i)=d*p(i)/r+x2(i) 1603 2 continue 1604c 1605 call rotate(x2,v,a,s,w) 1606 call rotate(x3,x2,t,w,s) 1607c 1608 do 3 i=1,3 1609 x1(i)=s(i) 1610 3 continue 1611c 1612 md_zmat=.true. 1613 return 1614 end 1615 subroutine md_abort(string, icode) 1616 implicit none 1617#include "global.fh" 1618#include "stdio.fh" 1619 character*(*) string 1620 character*255 card 1621 integer icode 1622 if(ga_nodeid().eq.0) then 1623 write(luout,1000) 0,string,icode 1624 1000 format(/,1x,10('*'),/,' * ',i3,': ',a,i5,/,1x,10('*')) 1625 card=' ' 1626 else 1627 write(card,1001) ga_nodeid(),string,icode 1628 1001 format(' * ',i3,': ',a,i5,' * ') 1629 endif 1630 call ga_error(card,icode) 1631 return 1632 end 1633