1module util 2implicit real*8(a-h,o-z) 3 4interface sort 5 module procedure sortr8 6 module procedure sorti4 7end interface 8interface invarr 9 module procedure invarrr8 10 module procedure invarri4 11end interface 12!!------------------- Root and weight of hermite polynomial 13real*8 Rhm(10,10),Whm(10,10) 14data Rhm(1,1) / 0.0D0 / 15data Rhm(2,1) / -0.70710678118654752440D+00 / 16data Rhm(2,2) / 0.70710678118654752440D+00 / 17data Rhm(3,1) / -1.22474487139158904910D+00 / 18data Rhm(3,2) / 0.0D0 / 19data Rhm(3,3) / 1.22474487139158904910D+00 / 20data Rhm(4,1) / -1.65068012388578455588D+00 / 21data Rhm(4,2) / -0.52464762327529031788D+00 / 22data Rhm(4,3) / 0.52464762327529031788D+00 / 23data Rhm(4,4) / 1.65068012388578455588D+00 / 24data Rhm(5,1) / -2.02018287045608563293D+00 / 25data Rhm(5,2) / -0.95857246461381850711D+00 / 26data Rhm(5,3) / 0.0D0 / 27data Rhm(5,4) / 0.95857246461381850711D+00 / 28data Rhm(5,5) / 2.02018287045608563293D+00 / 29data Rhm(6,1) / -2.35060497367449222283D+00 / 30data Rhm(6,2) / -1.33584907401369694971D+00 / 31data Rhm(6,3) / -0.43607741192761650868D+00 / 32data Rhm(6,4) / 0.43607741192761650868D+00 / 33data Rhm(6,5) / 1.33584907401369694971D+00 / 34data Rhm(6,6) / 2.35060497367449222283D+00 / 35data Rhm(7,1) / -2.65196135683523349245D+00 / 36data Rhm(7,2) / -1.67355162876747144503D+00 / 37data Rhm(7,3) / -0.81628788285896466304D+00 / 38data Rhm(7,4) / 0.0D0 / 39data Rhm(7,5) / 0.81628788285896466304D+00 / 40data Rhm(7,6) / 1.67355162876747144503D+00 / 41data Rhm(7,7) / 2.65196135683523349245D+00 / 42data Rhm(8,1) / -2.93063742025724401922D+00 / 43data Rhm(8,2) / -1.98165675669584292585D+00 / 44data Rhm(8,3) / -1.15719371244678019472D+00 / 45data Rhm(8,4) / -0.38118699020732211685D+00 / 46data Rhm(8,5) / 0.38118699020732211685D+00 / 47data Rhm(8,6) / 1.15719371244678019472D+00 / 48data Rhm(8,7) / 1.98165675669584292585D+00 / 49data Rhm(8,8) / 2.93063742025724401922D+00 / 50data Rhm(9,1) / -3.19099320178152760723D+00 / 51data Rhm(9,2) / -2.26658058453184311180D+00 / 52data Rhm(9,3) / -1.46855328921666793167D+00 / 53data Rhm(9,4) / -0.72355101875283757332D+00 / 54data Rhm(9,5) / 0.0D0 / 55data Rhm(9,6) / 0.72355101875283757332D+00 / 56data Rhm(9,7) / 1.46855328921666793167D+00 / 57data Rhm(9,8) / 2.26658058453184311180D+00 / 58data Rhm(9,9) / 3.19099320178152760723D+00 / 59data Rhm(10,1) / -3.43615911883773760333D+00 / 60data Rhm(10,2) / -2.53273167423278979641D+00 / 61data Rhm(10,3) / -1.75668364929988177345D+00 / 62data Rhm(10,4) / -1.03661082978951365418D+00 / 63data Rhm(10,5) / -0.34290132722370460879D+00 / 64data Rhm(10,6) / 0.34290132722370460879D+00 / 65data Rhm(10,7) / 1.03661082978951365418D+00 / 66data Rhm(10,8) / 1.75668364929988177345D+00 / 67data Rhm(10,9) / 2.53273167423278979641D+00 / 68data Rhm(10,10) / 3.43615911883773760333D+00 / 69data Whm(1,1) / 1.77245385090551602730D+00 / ! SQRT(PI) 70data Whm(2,1) / 8.86226925452758013649D-01 / 71data Whm(2,2) / 8.86226925452758013649D-01 / 72data Whm(3,1) / 2.95408975150919337883D-01 / 73data Whm(3,2) / 1.18163590060367735153D+00 / 74data Whm(3,3) / 2.95408975150919337883D-01 / 75data Whm(4,1) / 8.13128354472451771430D-02 / 76data Whm(4,2) / 8.04914090005512836506D-01 / 77data Whm(4,3) / 8.04914090005512836506D-01 / 78data Whm(4,4) / 8.13128354472451771430D-02 / 79data Whm(5,1) / 1.99532420590459132077D-02 / 80data Whm(5,2) / 3.93619323152241159828D-01 / 81data Whm(5,3) / 9.45308720482941881226D-01 / 82data Whm(5,4) / 3.93619323152241159828D-01 / 83data Whm(5,5) / 1.99532420590459132077D-02 / 84data Whm(6,1) / 4.53000990550884564086D-03 / 85data Whm(6,2) / 1.57067320322856643916D-01 / 86data Whm(6,3) / 7.24629595224392524092D-01 / 87data Whm(6,4) / 7.24629595224392524092D-01 / 88data Whm(6,5) / 1.57067320322856643916D-01 / 89data Whm(6,6) / 4.53000990550884564086D-03 / 90data Whm(7,1) / 9.71781245099519154149D-04 / 91data Whm(7,2) / 5.45155828191270305922D-02 / 92data Whm(7,3) / 4.25607252610127800520D-01 / 93data Whm(7,4) / 8.10264617556807326765D-01 / 94data Whm(7,5) / 4.25607252610127800520D-01 / 95data Whm(7,6) / 5.45155828191270305922D-02 / 96data Whm(7,7) / 9.71781245099519154149D-04 / 97data Whm(8,1) / 1.99604072211367619206D-04 / 98data Whm(8,2) / 1.70779830074134754562D-02 / 99data Whm(8,3) / 2.07802325814891879543D-01 / 100data Whm(8,4) / 6.61147012558241291030D-01 / 101data Whm(8,5) / 6.61147012558241291030D-01 / 102data Whm(8,6) / 2.07802325814891879543D-01 / 103data Whm(8,7) / 1.70779830074134754562D-02 / 104data Whm(8,8) / 1.99604072211367619206D-04 / 105data Whm(9,1) / 3.96069772632643819046D-05 / 106data Whm(9,2) / 4.94362427553694721722D-03 / 107data Whm(9,3) / 8.84745273943765732880D-02 / 108data Whm(9,4) / 4.32651559002555750200D-01 / 109data Whm(9,5) / 7.20235215606050957124D-01 / 110data Whm(9,6) / 4.32651559002555750200D-01 / 111data Whm(9,7) / 8.84745273943765732880D-02 / 112data Whm(9,8) / 4.94362427553694721722D-03 / 113data Whm(9,9) / 3.96069772632643819046D-05 / 114data Whm(10,1) / 7.64043285523262062916D-06 / 115data Whm(10,2) / 1.34364574678123269220D-03 / 116data Whm(10,3) / 3.38743944554810631362D-02 / 117data Whm(10,4) / 2.40138611082314686417D-01 / 118data Whm(10,5) / 6.10862633735325798784D-01 / 119data Whm(10,6) / 6.10862633735325798784D-01 / 120data Whm(10,7) / 2.40138611082314686417D-01 / 121data Whm(10,8) / 3.38743944554810631362D-02 / 122data Whm(10,9) / 1.34364574678123269220D-03 / 123data Whm(10,10) / 7.64043285523262062916D-06 / 124 125contains 126!Content sequences: 127!!Geometry operation 128!!Array, Vector 129!!String process 130!!Matrix calculation 131!!Misc 132 133!===============================================================! 134!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 135!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 136!!!!!!!!!!!!!!!!!!!!!!!!! Geometry operation !!!!!!!!!!!!!!!!!!!! 137!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 138!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 139!===============================================================! 140 141 142!!---------- Get angle (degree) between two vectors 143real*8 function vecang(vec1x,vec1y,vec1z,vec2x,vec2y,vec2z) 144real*8 vec1x,vec1y,vec1z,vec2x,vec2y,vec2z 145pi=3.141592653589793D0 146rnorm1=dsqrt(vec1x**2+vec1y**2+vec1z**2) 147rnorm2=dsqrt(vec2x**2+vec2y**2+vec2z**2) 148costheta=(vec1x*vec2x+vec1y*vec2y+vec1z*vec2z)/rnorm1/rnorm2 149if (costheta>1D0) costheta=1 150vecang=acos(costheta)/pi*180 151end function 152 153!!---------- Get distance of point 0 to plane(defined by point 1,2,3) 154real*8 function potpledis(x1,y1,z1,x2,y2,z2,x3,y3,z3,x0,y0,z0) 155real*8 x1,y1,z1,x2,y2,z2,x3,y3,z3,x0,y0,z0,prjx,prjy,prjz 156call pointprjple(x1,y1,z1,x2,y2,z2,x3,y3,z3,x0,y0,z0,prjx,prjy,prjz) 157potpledis=dsqrt((x0-prjx)**2+(y0-prjy)**2+(z0-prjz)**2) 158end function 159 160!!---------- Project a point(x0,y0,z0) to a plane defined by x/y/z-1/2/3, prjx/y/z are results 161subroutine pointprjple(x1,y1,z1,x2,y2,z2,x3,y3,z3,x0,y0,z0,prjx,prjy,prjz) 162real*8 x1,y1,z1,x2,y2,z2,x3,y3,z3,x0,y0,z0,prjx,prjy,prjz,A,B,C,D,t 163call pointABCD(x1,y1,z1,x2,y2,z2,x3,y3,z3,A,B,C,D) 164! (x0-x)/A=(y0-y)/B=(z0-z)/C ---> x=x0-t*A y=y0-t*B z=z0-t*C , substitute into Ax+By+Cz+D=0 solve t 165t=(D+A*x0+B*y0+C*z0)/(A**2+B**2+C**2) 166prjx=x0-t*A 167prjy=y0-t*B 168prjz=z0-t*C 169end subroutine 170 171!!-------- Input three points, get ABCD of Ax+By+Cz+D=0 172subroutine pointABCD(x1,y1,z1,x2,y2,z2,x3,y3,z3,A,B,C,D) 173real*8 v1x,v1y,v1z,v2x,v2y,v2z,x1,y1,z1,x2,y2,z2,x3,y3,z3,A,B,C,D 174v1x=x2-x1 175v1y=y2-y1 176v1z=z2-z1 177v2x=x3-x1 178v2y=y3-y1 179v2z=z3-z1 180! Solve determinant(Vector multiply) to get the normal vector (A,B,C): 181! i j k //unit vector 182! v1x v1y v1z 183! v2x v2y v2z 184A=v1y*v2z-v1z*v2y 185B=-(v1x*v2z-v1z*v2x) 186C=v1x*v2y-v1y*v2x 187D=A*(-x1)+B*(-y1)+C*(-z1) 188end subroutine 189 190!!-------- Input three points 0,1,2, get the vertical projection point of 0 to the line linking 1-2 191subroutine pointprjline(x0,y0,z0,x1,y1,z1,x2,y2,z2,prjx,prjy,prjz) 192real*8 x0,y0,z0,x1,y1,z1,x2,y2,z2,prjx,prjy,prjz,v12x,v12y,v12z,t 193v12x=x2-x1 194v12y=y2-y1 195v12z=z2-z1 196!Since prjx=x1+t*v12x prjy=y1+t*v12y prjz=z1+t*v12z 197!So v12x*(x0-prjx)+v12y*(y0-prjy)+v12z*(z0-prjz)=0 198!v12x*(x0-x1)-v12x*v12x*t + v12y*(y0-y1)-v12y*v12y*t + v12z*(z0-z1)-v12z*v12z*t =0 199t=( v12x*(x0-x1)+v12y*(y0-y1)+v12z*(z0-z1) )/(v12x**2+v12y**2+v12z**2) 200prjx=x1+t*v12x 201prjy=y1+t*v12y 202prjz=z1+t*v12z 203end subroutine 204 205!!---------- Get distance of point 0 to the line 1-2 206real*8 function potlinedis(x0,y0,z0,x1,y1,z1,x2,y2,z2) 207real*8 x0,y0,z0,x1,y1,z1,x2,y2,z2,prjx,prjy,prjz 208call pointprjline(x0,y0,z0,x1,y1,z1,x2,y2,z2,prjx,prjy,prjz) 209potlinedis=dsqrt((x0-prjx)**2+(y0-prjy)**2+(z0-prjz)**2) 210end function 211 212!!--------- Input three points, return angle between 1-2 and 2-3 (in degree) 213real*8 function xyz2angle(x1,y1,z1,x2,y2,z2,x3,y3,z3) 214real*8 :: x1,y1,z1,x2,y2,z2,x3,y3,z3,pi=3.141592653589793D0 215vec1x=x1-x2 216vec1y=y1-y2 217vec1z=z1-z2 218vec2x=x3-x2 219vec2y=y3-y2 220vec2z=z3-z2 221dotprod=vec1x*vec2x+vec1y*vec2y+vec1z*vec2z 222rnormv1=dsqrt( vec1x**2+vec1y**2+vec1z**2 ) 223rnormv2=dsqrt( vec2x**2+vec2y**2+vec2z**2 ) 224xyz2angle=acos(dotprod/(rnormv1*rnormv2))/pi*180 225end function 226 227!!--------- Input four points, return dihedral angle (in degree) 228real*8 function xyz2dih(x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4) 229real*8 :: x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,pi=3.141592653589793D0,phi 230v12x=x1-x2 231v12y=y1-y2 232v12z=z1-z2 233v23x=x2-x3 234v23y=y2-y3 235v23z=z2-z3 236v34x=x3-x4 237v34y=y3-y4 238v34z=z3-z4 239call vecprod(v12x,v12y,v12z,v23x,v23y,v23z,p1x,p1y,p1z) 240call vecprod(v23x,v23y,v23z,v34x,v34y,v34z,p2x,p2y,p2z) 241phi=acos( (p1x*p2x+p1y*p2y+p1z*p2z)/(sqrt(p1x*p1x+p1y*p1y+p1z*p1z)*sqrt(p2x*p2x+p2y*p2y+p2z*p2z)) ) 242xyz2dih=phi/pi*180 243end function 244 245!!--------------- Get area of a triangle, need input coordinates of three points 246function gettriangarea(pax,pay,paz,pbx,pby,pbz,pcx,pcy,pcz) 247implicit real*8 (a-h,o-z) 248real*8 gettriangarea,pax,pay,paz,pbx,pby,pbz,pcx,pcy,pcz 249! a---b va=pb-pa vb=pc-pa 250! | 251! V 252! c 253va1=pbx-pax 254va2=pby-pay 255va3=pbz-paz 256vb1=pcx-pax 257vb2=pcy-pay 258vb3=pcz-paz 259call vecprod(va1,va2,va3,vb1,vb2,vb3,vc1,vc2,vc3) !vc=va��vb=|va||vb|sin��*i where i is unit vector perpendicular to va and vb 260absvc=dsqrt(vc1**2+vc2**2+vc3**2) 261gettriangarea=0.5D0*absvc 262end function 263 264!!--------------- Get volume of a tetrahedron, need input coordinates of four points 265function gettetravol(pax,pay,paz,pbx,pby,pbz,pcx,pcy,pcz,pdx,pdy,pdz) 266implicit real*8 (a-h,o-z) 267real*8 pax,pay,paz,pbx,pby,pbz,pcx,pcy,pcz,pdx,pdy,pdz 268! real*8 volmat(4,4) 269! volmat(:,1)=(/ pax,pbx,pcx,pdx /) 270! volmat(:,2)=(/ pay,pby,pcy,pdy /) 271! volmat(:,3)=(/ paz,pbz,pcz,pdz /) 272! volmat(:,4)=1D0 273! gettetravol=abs(detmat(volmat))/6D0 274! call showmatgau(volmat) 275!vol=abs( (a-d)��((b-d)��(c-d)) )/6, see http://en.wikipedia.org/wiki/Tetrahedron 276vec1x=pax-pdx 277vec1y=pay-pdy 278vec1z=paz-pdz 279vec2x=pbx-pdx 280vec2y=pby-pdy 281vec2z=pbz-pdz 282vec3x=pcx-pdx 283vec3y=pcy-pdy 284vec3z=pcz-pdz 285call vecprod(vec2x,vec2y,vec2z,vec3x,vec3y,vec3z,vec2x3x,vec2x3y,vec2x3z) 286gettetravol=abs(vec1x*vec2x3x+vec1y*vec2x3y+vec1z*vec2x3z)/6D0 287end function 288 289 290 291 292!===============================================================! 293!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 294!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 295!!!!!!!!!!!!!!!!!!!!!!!! Array, Vector !!!!!!!!!!!!!!!!!!!!!!!!!! 296!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 297!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 298!===============================================================! 299 300!!--------- Invert integer array, from istart to iend 301!Real*8 version 302subroutine invarrr8(array,istart,iend) 303integer istart,iend 304real*8 array(:) 305len=iend-istart+1 306do i=0,int(len/2D0)-1 307 tmp=array(istart+i) 308 array(istart+i)=array(iend-i) 309 array(iend-i)=tmp 310end do 311end subroutine 312!Integer 4 version 313subroutine invarri4(array,istart,iend) 314integer istart,iend,array(:) 315len=iend-istart+1 316do i=0,int(len/2D0)-1 317 tmp=array(istart+i) 318 array(istart+i)=array(iend-i) 319 array(iend-i)=tmp 320end do 321end subroutine 322 323!-------- Sort value from small to big by Bubble method 324!mode=abs: sort by absolute value, =val: sort by value. Default is by value 325!Real*8 version 326subroutine sortr8(array,inmode) 327integer N,i,j,mode 328character,optional :: inmode*3 329real*8 array(:),temp 330N=size(array) 331mode=1 332if (present(inmode)) then 333 if (inmode=="abs") mode=2 334end if 335if (mode==1) then 336 do i=1,N 337 do j=i+1,N 338 if (array(i)>array(j)) then 339 temp=array(i) 340 array(i)=array(j) 341 array(j)=temp 342 end if 343 end do 344 end do 345else if (mode==2) then 346 do i=1,N 347 do j=i+1,N 348 if (abs(array(i))>abs(array(j))) then 349 temp=array(i) 350 array(i)=array(j) 351 array(j)=temp 352 end if 353 end do 354 end do 355end if 356end subroutine 357!Integer 4 version 358subroutine sorti4(array,inmode) 359integer N,i,j,mode 360character,optional :: inmode*3 361integer array(:),temp 362N=size(array) 363mode=1 364if (present(inmode)) then 365 if (inmode=="abs") mode=2 366end if 367if (mode==1) then 368 do i=1,N 369 do j=i+1,N 370 if (array(i)>array(j)) then 371 temp=array(i) 372 array(i)=array(j) 373 array(j)=temp 374 end if 375 end do 376 end do 377else if (mode==2) then 378 do i=1,N 379 do j=i+1,N 380 if (abs(array(i))>abs(array(j))) then 381 temp=array(i) 382 array(i)=array(j) 383 array(j)=temp 384 end if 385 end do 386 end do 387end if 388end subroutine 389 390!!--------- Evaluate standard deviation of array elements 391real*8 function stddevarray(array) 392real*8 array(:),avg 393avg=sum(array)/size(array) 394stddevarray=dsqrt(sum((array-avg)**2)/size(array)) 395end function 396 397!!--------- Evaluate covariant of two array elements 398real*8 function covarray(array1,array2) 399real*8 array1(:),array2(:),avg1,avg2 400avg1=sum(array1)/size(array1) 401avg2=sum(array2)/size(array2) 402covarray=sum((array1-avg1)*(array2-avg2))/size(array1) 403end function 404 405!--- Vector/cross product, input two vectors, return a new vector (x,y,z) 406subroutine vecprod(x1,y1,z1,x2,y2,z2,x,y,z) 407real*8 x1,y1,z1,x2,y2,z2,x,y,z 408! |i j k | 409! |x1 y1 z1| 410! |x2 y2 z2| 411x= y1*z2-z1*y2 412y=-(x1*z2-z1*x2) 413z= x1*y2-y1*x2 414end subroutine 415 416!--- Generate full arrangement arry 417!ncol is the number of element, nrow=ncol! 418!example: nrow=3!=5, ncol=3, the arr will be: 419! 1 3 2 420! 2 1 3 421! 2 3 1 422! 3 1 2 423! 3 2 1 424! 1 2 3 425subroutine fullarrange(arr,nrow,ncol) 426integer nrow,ncol,arr(nrow,ncol),seq(ncol) 427seq=(/ (i,i=1,ncol) /) 428arr(1,:)=seq !The first array will be 1,2,3,4... 429do icyc=2,nrow 430 do i=ncol-1,1,-1 431 if (seq(i)<seq(i+1)) exit 432 end do 433 do j=ncol,1,-1 434 if (seq(j)>seq(i)) exit 435 end do 436 itmp=seq(i) 437 seq(i)=seq(j) 438 seq(j)=itmp 439 call invarr(seq,i+1,ncol) 440 arr(icyc,:)=seq 441end do 442end subroutine 443 444 445 446!===============================================================! 447!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 448!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 449!!!!!!!!!!!!!!!!!!!!!!!!!!! String process !!!!!!!!!!!!!!!!!!!!!! 450!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 451!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 452!===============================================================! 453 454!-------- Parse inputted integer string (e.g. 3,4,5,6,7-25,32-35,55,88) to array 455!nelement will be return, which is the number of indices in the string 456!if "array" is present, the terms will be written into it, else if "array" is not present, only count the number of terms as "nelement" 457!array must be large enough to contain the elements 458subroutine str2arr(inpstr,nelement,array) 459character(len=*) inpstr 460character c80tmp*80 461integer nelement 462integer,optional :: array(:) 463if (present(array)) array=0 464icommaold=0 465nelement=0 466do ipos=1,len(inpstr) 467 if (inpstr(ipos:ipos)==','.or.ipos==len(inpstr)) then 468 icommanew=ipos 469 if (ipos==len(inpstr)) icommanew=ipos+1 470 read(inpstr(icommaold+1:icommanew-1),*) c80tmp 471 if (index(c80tmp,'-')==0) then 472 nelement=nelement+1 473 if (present(array)) read(c80tmp,*) array(nelement) 474 else 475 do iheng=1,len_trim(c80tmp) 476 if (c80tmp(iheng:iheng)=='-') exit 477 end do 478 read(c80tmp(1:iheng-1),*) ilow 479 read(c80tmp(iheng+1:),*) ihigh 480 do itmp=ilow,ihigh 481 nelement=nelement+1 482 if (present(array)) array(nelement)=itmp 483 end do 484 end if 485 icommaold=icommanew 486 end if 487end do 488end subroutine 489 490!---------Input path name, e.g. c:\ltwd\MIO.wfn , output file name, e.g. MIO 491subroutine path2filename(pathnamein,filenameout) 492character(len=*) pathnamein,filenameout 493do i=len_trim(pathnamein),1,-1 494 if (pathnamein(i:i)=='.') then 495 iend=i-1 496 exit 497 end if 498end do 499istart=1 500do i=iend,1,-1 501 if (pathnamein(i:i)=='/'.or.pathnamein(i:i)=='\') then 502 istart=i+1 503 exit 504 end if 505end do 506filenameout=' ' 507filenameout(1:iend-istart+1)=pathnamein(istart:iend) 508end subroutine 509 510!!--------- Convert a character to lower case 511subroutine uc2lc(inc) 512character*1 inc 513itmp=ichar(inc) 514if (itmp>=65.and.itmp<=90) itmp=itmp+32 515inc=char(itmp) 516end subroutine 517 518!!--------- Convert a character to upper case 519subroutine lc2uc(inc) 520character*1 inc 521itmp=ichar(inc) 522if (itmp>=97.and.itmp<=122) itmp=itmp-32 523inc=char(itmp) 524end subroutine 525 526!!--------- Convert a string to lower case 527subroutine struc2lc(str) 528character(len=*) str 529do i=1,len_trim(str) 530 call uc2lc(str(i:i)) 531end do 532end subroutine 533 534!!--------- Convert a string to upper case 535subroutine strlc2uc(str) 536character(len=*) str 537do i=1,len_trim(str) 538 call lc2uc(str(i:i)) 539end do 540end subroutine 541 542 543 544!===============================================================! 545!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 546!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 547!!!!!!!!!!!!!!!!!!!!!!!!! Matrix calculation !!!!!!!!!!!!!!!!!!!! 548!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 549!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 550!===============================================================! 551 552 553!!----- Make the matrix to upper trigonal matrix 554subroutine ratio_upper(mat) 555real*8 :: mat(:,:),m,st 556real*8,allocatable :: temp(:),s(:),divided(:) 557integer :: a,i,j,t 558a=size(mat,1) 559allocate(temp(a)) 560allocate(s(a)) 561allocate(divided(a)) 562do i=1,a 563 s(i)=maxval(abs(mat(i,1:a))) 564end do 565do i=1,a-1 566 divided(i:a)=mat(i:a,i)/s(i:a) 567 t=maxloc(abs(divided(i:a)),dim=1) 568 temp(:)=mat(i,:) 569 mat(i,:)=mat(i+t-1,:) 570 mat(i+t-1,:)=temp(:) 571 st=s(i) 572 s(i)=s(i+t-1) 573 s(i+t-1)=st 574 do j=i+1,a 575 m=mat(j,i)/mat(i,i) 576 mat(j,i:a)=mat(j,i:a)-mat(i,i:a)*m 577 end do 578end do 579deallocate(temp,s,divided) 580end subroutine 581 582!!----- Get value of determinant of a matrix 583real*8 function detmat(mat) 584real*8 mat(:,:) 585real*8,allocatable :: mattmp(:,:) 586isizemat=size(mat,1) 587detmat=1D0 588NOTlowertri=0 589NOTuppertri=0 590outter1: do i=1,isizemat !Check if already is lower-trigonal matrix 591 do j=i+1,isizemat 592 if (mat(i,j)>1D-12) then 593 NOTlowertri=1 !There are at least one big value at upper trigonal part, hence not lower trigonal matrix 594 exit outter1 595 end if 596 end do 597end do outter1 598outter2: do i=1,isizemat !Check if already is upper-trigonal matrix 599 do j=1,i-1 600 if (mat(i,j)>1D-12) then 601 NOTuppertri=1 !There are at least one big value at lower trigonal part, hence not upper trigonal matrix 602 exit outter2 603 end if 604 end do 605end do outter2 606 607if (NOTlowertri==0.or.NOTuppertri==0) then !Is lower or upper trigonal matrix, don't need to convert to trigonal matrix 608 do i=1,isizemat 609 detmat=detmat*mat(i,i) 610 end do 611else !Not upper or lower trigonal matrix 612 allocate(mattmp(isizemat,isizemat)) 613 mattmp=mat 614 call ratio_upper(mattmp) 615 detmat=1D0 616 do i=1,isizemat 617 detmat=detmat*mattmp(i,i) 618 end do 619end if 620end function 621 622!!-------- Get trace of a matrix 623real*8 function mattrace(mat) 624real*8 mat(:,:) 625mattrace=0 626do i=1,size(mat,1) 627 mattrace=mattrace+mat(i,i) 628end do 629end function 630 631!!--- Use Jacobi method to diagonalize matrix, simple, but much slower than diagsymat and diaggemat if the matrix is large 632subroutine diagmat(mat,S,eigval,inmaxcyc,inthres) 633! mat: input and will be diagonalized matrix, S:eigenvector matrix(columns correspond to vectors), eigval:eigenvalue vector 634! inmaxcyc: max cycle, inthres: expected threshold 635implicit real*8 (a-h,o-z) 636integer,optional :: inmaxcyc 637real*8,optional :: inthres 638real*8 thres,mat(:,:),S(:,:),eigval(:) 639real*8,allocatable :: R(:,:) 640n=size(mat,1) 641allocate(R(n,n)) 642maxcyc=200 643thres=1D-9 644if (present(inmaxcyc)) maxcyc=inmaxcyc 645if (present(inthres)) thres=inthres 646S=0 647do i=1,n 648 S(i,i)=1.0D0 649end do 650do k=1,maxcyc+1 651 R=0 652 do i=1,n 653 R(i,i)=1.0D0 654 end do 655 i=1 656 j=2 657 do ii=1,n 658 do jj=ii+1,n 659 if (abs(mat(ii,jj))>abs(mat(i,j))) then 660 i=ii 661 j=jj 662 end if 663 end do 664 end do 665 if (abs(mat(i,j))<thres) exit 666 if (k==maxcyc+1) write(*,*) "Matrix diagonalization exceed max cycle before converge" 667 phi=atan(2*mat(i,j)/(mat(i,i)-mat(j,j)))/2.0D0 668 R(i,i)=cos(phi) 669 R(j,j)=R(i,i) 670 R(i,j)=-sin(phi) 671 R(j,i)=-R(i,j) 672 mat=matmul(matmul(transpose(R),mat),R) 673 S=matmul(S,R) 674end do 675do i=1,n 676 eigval(i)=mat(i,i) 677end do 678end subroutine 679 680!!------------ Diagonalize a symmetry matrix 681!Repack the extremely complex "DSYEV" routine in lapack to terse form 682!if istat/=0, means error occurs 683subroutine diagsymat(mat,eigvecmat,eigvalarr,istat) 684integer istat 685real*8 mat(:,:),eigvecmat(:,:),eigvalarr(:) 686real*8,allocatable :: lworkvec(:) 687isize=size(mat,1) 688allocate(lworkvec(3*isize-1)) 689call DSYEV('V','U',isize,mat,isize,eigvalarr,lworkvec,3*isize-1,istat) 690eigvecmat=mat 691mat=0D0 692forall (i=1:isize) mat(i,i)=eigvalarr(i) 693end subroutine 694 695!!------------ Diagonalize a general matrix 696!Repack the extremal complex "DGEEV" routine in lapack to terse form 697!eigvecmat is right eigenvector matrix 698!eigvalarr is real part of eigenvalue, imaginary parts are discarded 699!if istat/=0, means error appears 700subroutine diaggemat(mat,eigvecmat,eigvalarr,istat) 701integer istat,lwork 702real*8 mat(:,:),eigvecmat(:,:),eigvalarr(:),tmpmat(1,1) 703real*8,allocatable :: lworkvec(:),eigvalimgarr(:) 704isize=size(mat,1) 705lwork=8*isize !4*isize is enough, but for better performance we use larger size 706allocate(lworkvec(lwork),eigvalimgarr(isize)) 707call DGEEV('N','V',isize,mat,isize,eigvalarr,eigvalimgarr,tmpmat,1,eigvecmat,isize,lworkvec,lwork,istat) 708mat=0D0 709forall (i=1:isize) mat(i,i)=eigvalarr(i) 710end subroutine 711 712!!--------------- A function to output inverted matrix, inputted matrix will not be affected. Essentially is a warpper of KROUT 713function invmat(mat,N) 714integer N,ierr 715real*8 :: mat(N,N),invmat(N,N),tmpvec(N) 716invmat=mat 717call KROUT(0,N,0,invmat,N,tmpvec,N,ierr) 718end function 719!!--------------- A routine to invert matrix. The inputted matrix will be taken placed by inverted matrix. Essentially is a warpper of KROUT 720subroutine invmatsub(mat,N) 721integer N,ierr 722real*8 :: mat(N,N),tmpvec(N) 723call KROUT(0,N,0,mat,N,tmpvec,N,ierr) 724end subroutine 725!Taken and adapted from krout.f, which can be downloaded at http://jblevins.org/mirror/amiller/ 726!----------------------------------------------------------------------- 727! CROUT PROCEDURE FOR INVERTING MATRICES AND SOLVING EQUATIONS 728!----------------------------------------------------------------------- 729! A IS A MATRIX OF ORDER N WHERE N IS GREATER THAN OR EQUAL TO 1. 730! IF MO = 0 THEN THE INVERSE OF A IS COMPUTED AND STORED IN A. 731! IF MO IS NOT 0 THEN THE INVERSE IS NOT COMPUTED. 732 733! IF M IS GREATER THAN 0 THEN B IS A MATRIX HAVING N ROWS AND M COLUMNS. 734! IN THIS CASE AX = B IS SOLVED AND THE SOLUTION X IS STORED IN B. 735! IF M=0 THEN THERE ARE NO EQUATIONS TO BE SOLVED. 736! N.B. B is passed as a VECTOR not as a matrix. 737 738! KA = THE LENGTH OF THE COLUMNS OF THE ARRAY A 739! KB = THE LENGTH OF THE COLUMNS OF THE ARRAY B (IF M > 0) 740 741! IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS. WHEN 742! THE ROUTINE TERMINATES IERR HAS ONE OF THE FOLLOWING VALUES ... 743! IERR = 0 THE REQUESTED TASK WAS PERFORMED. 744! IERR = -1 EITHER N, KA, OR KB IS INCORRECT. 745! IERR = K THE K-TH PIVOT ELEMENT IS 0. 746!----------------------------------------------------------------------- 747! Adapted from the routine KROUT in the NSWC Math. Library by Alan Miller 748! Latest revision - 3 August 1998 749subroutine KROUT(MO, N, M, A, KA, B, KB, IERR) 750implicit none 751integer, intent(in) :: MO 752integer, intent(in) :: N 753integer, intent(in) :: M 754real*8, intent(in out), dimension(:,:) :: A ! a(ka,n) 755integer, intent(in) :: KA 756real*8, intent(in out), dimension(:) :: B 757integer, intent(in) :: KB 758integer, intent(out) :: IERR 759integer, allocatable, dimension(:) :: INDX 760real*8, allocatable, dimension(:) :: TEMP 761integer :: I, J, JP1, K, KJ, KM1, KP1, L, LJ, MAXB, NJ, NMJ, NMK, NM1, ONEJ 762real*8 :: D, DSUM, P, T 763real*8, parameter :: ZERO = 0.0D0, ONE = 1.0D0 764if (N < 1 .or. KA < N) then 765 IERR = -1 766 return 767end if 768if (M > 0 .and. KB < N) then 769 IERR = -1 770 return 771end if 772IERR = 0 773if (N < 2) then 774 D = A(1,1) 775 if (D == ZERO) then 776 IERR = N 777 return 778 end if 779 if (MO == 0) A(1,1) = ONE / D 780 if (M <= 0) return 781 MAXB = KB*M 782 do KJ = 1,MAXB,KB 783 B(KJ) = B(KJ)/D 784 end do 785 return 786end if 787if (MO == 0) then 788 allocate( INDX(N-1), TEMP(N) ) 789end if 790NM1 = N - 1 791do K = 1,NM1 792 KP1 = K + 1 793 P = abs(A(K,K)) 794 L = K 795 do I = KP1,N 796 T = abs(A(I,K)) 797 if (P >= T) cycle 798 P = T 799 L = I 800 end do 801 if (P == ZERO) then 802 IERR = K 803 return 804 end if 805 P = A(L,K) 806 if (MO == 0) then 807 INDX(K) = L 808 end if 809 if (K /= L) then 810 do J = 1,N 811 T = A(K,J) 812 A(K,J) = A(L,J) 813 A(L,J) = T 814 end do 815 if (M > 0) then 816 KJ = K 817 LJ = L 818 do J = 1,M 819 T = B(KJ) 820 B(KJ) = B(LJ) 821 B(LJ) = T 822 KJ = KJ + KB 823 LJ = LJ + KB 824 end do 825 end if 826 end if 827 if (K <= 1) then 828 do J = KP1,N 829 A(K,J) = A(K,J)/P 830 end do 831 else 832 do J = KP1,N 833 DSUM = A(K,J) - dot_product( A(K,1:KM1), A(1:KM1,J) ) 834 A(K,J) = DSUM / P 835 end do 836 end if 837 do I = KP1,N 838 DSUM = A(I,KP1) - dot_product( A(I,1:K), A(1:K,KP1) ) 839 A(I,KP1) = DSUM 840 end do 841 KM1 = K 842end do 843if (A(N,N) == ZERO) then 844 IERR = N 845 return 846end if 847if (M > 0) then 848 MAXB = KB*M 849 do ONEJ = 1,MAXB,KB 850 KJ = ONEJ 851 B(KJ) = B(KJ)/A(1,1) 852 do K = 2,N 853 KJ = KJ + 1 854 DSUM = B(KJ) 855 KM1 = K - 1 856 LJ = ONEJ 857 do L = 1,KM1 858 DSUM = DSUM - A(K,L)*B(LJ) 859 LJ = LJ + 1 860 end do 861 B(KJ) = DSUM / A(K,K) 862 end do 863 end do 864 do NJ = N,MAXB,KB 865 KJ = NJ 866 do NMK = 1,NM1 867 K = N - NMK 868 LJ = KJ 869 KJ = KJ - 1 870 DSUM = B(KJ) 871 KP1 = K + 1 872 do L = KP1,N 873 DSUM = DSUM - A(K,L)*B(LJ) 874 LJ = LJ + 1 875 end do 876 B(KJ) = DSUM 877 end do 878 end do 879end if 880if (MO /= 0) return 881do J = 1,NM1 882 A(J,J) = ONE / A(J,J) 883 JP1 = J + 1 884 do I = JP1,N 885 DSUM = dot_product( A(I,J:I-1), A(J:I-1,J) ) 886 A(I,J) = -DSUM / A(I,I) 887 end do 888end do 889A(N,N) = ONE / A(N,N) 890do NMK = 1,NM1 891 K = N - NMK 892 KP1 = K + 1 893 do J = KP1,N 894 TEMP(J) = A(K,J) 895 A(K,J) = ZERO 896 end do 897 do J = 1,N 898 DSUM = A(K,J) - dot_product( TEMP(KP1:N), A(KP1:N,J) ) 899 A(K,J) = DSUM 900 end do 901end do 902do NMJ = 1,NM1 903 J = N - NMJ 904 K = INDX(J) 905 if (J == K) cycle 906 do I = 1,N 907 T = A(I,J) 908 A(I,J) = A(I,K) 909 A(I,K) = T 910 end do 911end do 912if (MO == 0) deallocate( INDX, TEMP ) 913end subroutine 914 915 916!------- Calculate how much is a matrix deviates from identity matrix 917!error=��[i,j]abs( abs(mat(i,j))-��(i,j) ) 918real*8 function identmaterr(mat) 919implicit real*8 (a-h,o-z) 920real*8 mat(:,:) 921nsize=size(mat,1) 922identmaterr=0D0 923do i=1,nsize 924 do j=1,nsize 925 if (i==j) then 926 identmaterr=identmaterr+abs(abs(mat(i,j))-1D0) 927 else 928 identmaterr=identmaterr+abs(mat(i,j)) 929 end if 930 end do 931end do 932end function 933 934 935!----- Convert a square matrix to an array. imode=1/2/3: Full matrix; Lower half matrix; Upper half matrix 936!For mode=1,2, "arr" should be nsize*(nsize+1)/2 937subroutine mat2arr(mat,arr,imode) 938implicit real*8 (a-h,o-z) 939real*8 mat(:,:),arr(:) 940nsize=size(mat,1) 941itmp=0 942if (imode==1) then !Full matrix 943 do i=1,nsize 944 do j=1,nsize 945 itmp=itmp+1 946 arr(itmp)=mat(i,j) 947 end do 948 end do 949else if (imode==2) then !Lower half matrix 950 do i=1,nsize 951 do j=1,i 952 itmp=itmp+1 953 arr(itmp)=mat(i,j) 954 end do 955 end do 956else !Upper half matrix 957 do i=1,nsize 958 do j=i,nsize 959 itmp=itmp+1 960 arr(itmp)=mat(i,j) 961 end do 962 end do 963end if 964end subroutine 965 966 967!===============================================================! 968!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 969!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 970!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Misc !!!!!!!!!!!!!!!!!!!!!!!!!!!!! 971!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 972!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 973!===============================================================! 974 975!!---------- Get current time in second, the difference between two times of invoking this routine is consumed wall clock time 976subroutine walltime(inow) 977character nowdate*20,nowtime*20 978integer inow 979call date_and_time(nowdate,nowtime) 980read(nowtime(1:2),*) inowhour 981read(nowtime(3:4),*) inowminute 982read(nowtime(5:6),*) inowsecond 983inow=inowhour*3600+inowminute*60+inowsecond 984end subroutine 985 986!!----- Find the position of specific value in cube 987subroutine findvalincub(cubfile,value,i,j,k) 988real*8 cubfile(:,:,:),value 989integer i,j,k 990do ii=1,size(cubfile,1) 991 do jj=1,size(cubfile,2) 992 do kk=1,size(cubfile,3) 993 if (cubfile(ii,jj,kk)==value) then 994 i=ii 995 j=jj 996 k=kk 997 end if 998 end do 999 end do 1000end do 1001end subroutine 1002 1003!!------ Display matrix similar to gaussian program 1004subroutine showmatgau(mat,label,insemi,form,fileid,useri1,useri2,inncol,titlechar) 1005!The number of columns is always 5 (unadjustable) 1006!"Label" is the title, if the content is empty, title will not be printed 1007!If semi==1, only lower and diagonal element will be shown 1008!"form" is the format to show data, default is D14.6, can pass into such as "f14.8", total width should be 14 characters 1009!fildid is output destination, 6 corresponds to outputting to screen 1010!useri1 and useri2 define the dimension of the matrix, default or -1 means determining automatically 1011!inncol seems controls spacing between number labels of each frame 1012!Default titlechar is "i8,6x", if you manually set inncol, you also set this to broaden or narrow 1013implicit real*8(a-h,o-z) 1014real*8 :: mat(:,:) 1015character(*),optional :: label,form,titlechar 1016integer,optional :: insemi,fileid,useri1,useri2,inncol 1017integer :: semi,ides,ncol 1018semi=0 1019ides=6 1020ncol=5 1021i1=size(mat,1) 1022i2=size(mat,2) 1023if (present(useri1).and.useri1/=-1) i1=useri1 1024if (present(useri2).and.useri1/=-1) i2=useri2 1025if (present(insemi)) semi=insemi 1026if (present(fileid)) ides=fileid 1027if (present(inncol)) ncol=inncol 1028if (present(label).and.label/='') write(ides,*) "************ ",label," ************" 1029nt=ceiling(i2/float(ncol)) 1030do i=1,nt !How many frame 1031 ns=(i-1)*5+1 !This frame start from where 1032 if (i/=nt) ne=(i-1)*ncol+ncol !This frame end to where 1033 if (i==nt) ne=i2 1034 !Write basis number in separate line 1035 write(ides,"(6x)",advance='no') 1036 do j=ns,ne 1037 if (present(titlechar)) then 1038 write(ides,'('//titlechar//')',advance='no') j 1039 else 1040 write(ides,"(i8,6x)",advance='no') j 1041 end if 1042 end do 1043 write(ides,*) 1044 !Write content in each regular line 1045 do k=1,i1 1046 if (k<ns.and.semi==1) cycle !The lines have been outputted are skipped 1047 write(ides,"(i6)",advance='no') k 1048 do j=ns,ne 1049 if (semi==1.and.k<j) cycle !Upper trigonal element were passed 1050 if (present(form)) then 1051 write(ides,'('//form//')',advance='no') mat(k,j) 1052 else 1053 write(ides,"(D14.6)",advance='no') mat(k,j) 1054 end if 1055 end do 1056 write(ides,*) !Change to next line 1057 end do 1058end do 1059end subroutine 1060 1061!!------- Read matrix outputted in the gaussian .out file, such as outputted by iop(3/33=1) 1062subroutine readmatgau(fileid,mat,semi,inform,inskipcol,inncol,innspace) 1063!e.g. trigonal: readmatgau(10,tempmat,1,"D14.6",7,5) full: readmatgau(10,tempmat,0,"D13.6",7,5) 1064!inform is the format used to read data, default is D14.6, you can input "f8.4 " (Note: Must be 5 character!) 1065!semi=1 means the read matrix is lower trigonal matrix 1066!inskipcol is the skipped column (marker information) in front of each row, default is 7 1067!inncol is data in each row, default is 5 1068!innspace is space lines in between each frame, default is 1 1069 1070!Before use, loclabel should be used to move reading position into title line of the matrix��namely move to "*** Overlap ***" 1071!This routine will skip title line, and skip one line in each frame reading 1072! *** Overlap *** 1073! 1 2 3 4 5 1074! 1 0.100000D+01 1075! 2 0.236704D+00 0.100000D+01 1076implicit real*8(a-h,o-z) 1077real*8 :: mat(:,:) 1078character,optional :: inform*5 1079character :: form*7,c80tmp*79 1080integer,optional :: inskipcol,inncol,semi,innspace 1081integer :: skipcol,ncol,fileid 1082form="(D14.6)" !Suitable for Sbas,Kinene,Potene,Hcore by 3/33=1 ,Fockmat,Densmat by 5/33=3 1083skipcol=7 1084ncol=5 1085nspace=1 1086if (present(inform)) form='('//inform//')' 1087if (present(inskipcol)) skipcol=inskipcol 1088if (present(inncol)) ncol=inncol 1089if (present(innspace)) nspace=innspace 1090read(fileid,*) !!!!!!!!!!!! Skip title line 1091i1=size(mat,1) 1092i2=size(mat,2) 1093nt=ceiling(i2/float(ncol)) 1094mat=0D0 1095do i=1,nt !Number of frames 1096 ns=(i-1)*ncol+1 1097 if (i/=nt) ne=(i-1)*ncol+ncol 1098 if (i==nt) ne=i2 1099 do ii=1,nspace 1100 read(fileid,*) !!!!!!!!!!!! Skip number line when reading each new frame 1101 end do 1102 do k=1,i1 !Scan rows in each frame 1103! read(fileid,"(a)") c80tmp 1104! write(15,"(a)") c80tmp 1105! backspace(fileid) 1106 1107 if (k<ns.and.present(semi).and.semi==1) cycle 1108 do iii=1,skipcol !Skip marker columns in each row 1109 read(fileid,"(1x)",advance='no') 1110 end do 1111 do j=ns,ne !Scan elements in each row 1112 if (present(semi).and.semi==1.and.k<j) cycle 1113 read(fileid,form,advance='no') mat(k,j) 1114! write(*,*) i,k,j,mat(k,j) 1115 end do 1116 read(fileid,*) 1117 end do 1118end do 1119if (present(semi).and.semi==1) then !When read is lower trigonal matrix, we assume it is a symmetric matrix 1120 mat=mat+transpose(mat) 1121 do i=1,i1 1122 mat(i,i)=mat(i,i)/2D0 1123 end do 1124end if 1125end subroutine 1126 1127!!!------------------------- Determine how many lines in the fileid 1128!If imode==1, space line will be regarded as the sign of end of file. If imode==2, will count actual number of lines in the file 1129integer function totlinenum(fileid,imode) 1130integer fileid,ierror,imode 1131character*80 c80 1132totlinenum=0 1133do while(.true.) 1134 read(fileid,"(a)",iostat=ierror) c80 1135 if (imode==1) then 1136 if (ierror/=0.or.c80==" ") exit 1137 else if (imode==2) then 1138 if (ierror/=0) exit 1139 end if 1140 totlinenum=totlinenum+1 1141end do 1142rewind(fileid) 1143end function 1144 1145!!!-------- Locate the line where the label first appears in fileid 1146!Return ifound=1 if found the label, else return 0 1147!Default is rewind, if irewind=0 then will not rewind 1148!If the current line just has the label, calling this subroutine will do nothing 1149!maxline define the maximum number of lines that will be searched, default is search the whole file 1150subroutine loclabel(fileid,label,ifound,irewind,maxline) 1151integer fileid,ierror 1152integer,optional :: ifound,irewind,maxline 1153character*200 c200 1154CHARACTER(LEN=*) label 1155if ((.not.present(irewind)).or.(present(irewind).and.irewind==1)) rewind(fileid) 1156if (.not.present(maxline)) then 1157 do while(.true.) 1158 read(fileid,"(a)",iostat=ierror) c200 1159 if (index(c200,label)/=0) then 1160 backspace(fileid) 1161 if (present(ifound)) ifound=1 !Found result 1162 return 1163 end if 1164 if (ierror/=0) exit 1165 end do 1166else 1167 do iline=1,maxline 1168 read(fileid,"(a)",iostat=ierror) c200 1169 if (index(c200,label)/=0) then 1170 backspace(fileid) 1171 if (present(ifound)) ifound=1 !Found result 1172 return 1173 end if 1174 if (ierror/=0) exit 1175 end do 1176end if 1177if (present(ifound)) ifound=0 1178end subroutine 1179 1180 1181!!!----------- Skip specific number of lines in specific fileid 1182subroutine skiplines(id,nskip) 1183integer id,nskip 1184do i=1,nskip 1185 read(id,*) 1186end do 1187end subroutine 1188 1189 1190!!!---------------- Calculate factorial 1191integer function ft(i) 1192integer i 1193ft=i 1194if (i==0) ft=1 1195do j=i-1,1,-1 1196 ft=ft*j 1197end do 1198end function 1199 1200!---- Calculate gamma(Lval+1/2), see http://en.wikipedia.org/wiki/Gamma_function 1201real*8 function gamma_ps(n) 1202use defvar 1203integer n 1204gamma_ps=ft(2*n)*dsqrt(pi)/4**n/ft(n) 1205end function 1206 1207!!-------- Get all combinations of any ncomb elements of array, which length is ntot 1208!outarray(A,B) is output array, A is the length and must be ntot!/ncomb!/(ntot-ncomb)!, B is generated array, should be equal to ncomb 1209subroutine combarray(array,ntot,ncomb,outarray) 1210integer array(ntot),ntot,ncomb,idxarr(ncomb),outarray(:,:) 1211ipos=ncomb !Current position in the array 1212forall (i=1:ncomb) idxarr(i)=i !Used to record index 1213ioutput=1 1214ncount=0 1215do while(ipos>0) 1216 if (ioutput==1) then 1217 ncount=ncount+1 1218 outarray(ncount,:)=array(idxarr(:)) 1219 end if 1220 ioutput=0 1221 idxarr(ipos)=idxarr(ipos)+1 1222 if (idxarr(ipos)>ntot) then 1223 ipos=ipos-1 !Go back to last position 1224 cycle 1225 end if 1226 if (ipos<ncomb) then 1227 ipos=ipos+1 1228 idxarr(ipos)=idxarr(ipos-1) 1229 cycle 1230 end if 1231 if (ipos==ncomb) ioutput=1 1232end do 1233end subroutine 1234 1235 1236!!--------- Convert XY scatter data to density distribution 1237!xarr and yarr records the points. nlen is array length 1238!mat is the outputted matrix, matnx and matny are its number of element in X and Y 1239!x/ymin, x/ymax are lower and upper limit of "mat", the X and Y ranges contain nvalx and nvaly data 1240!e.g. n=5 1241! | 1 | 2 | 3 | 4 | 5 | 1242! xmin-------------------------------------xmax 1243subroutine xypt2densmat(xarr,yarr,nlen,mat,nvalx,nvaly,xmin,xmax,ymin,ymax) 1244integer nlen,nvalx,nvaly 1245real*8 xarr(nlen),yarr(nlen),mat(nvalx,nvaly),xmin,xmax,ymin,ymax 1246mat=0D0 1247spcx=(xmax-xmin)/nvalx 1248spcy=(ymax-ymin)/nvaly 1249!If enable parallel, program often prompts memory is not enough, I don't know how to solve this 1250! !$OMP PARALLEL DO SHARED(mat) PRIVATE(ix,iy,ipt,xlow,xhigh,ylow,yhigh) schedule(dynamic) NUM_THREADS(nthreads) 1251do ix=1,nvalx 1252 xlow=xmin+(ix-1)*spcx 1253 xhigh=xmin+ix*spcx 1254 do iy=1,nvaly 1255 ylow=ymin+(iy-1)*spcy 1256 yhigh=ymin+iy*spcy 1257 do ipt=1,nlen 1258 if (xarr(ipt)>xlow.and.xarr(ipt)<=xhigh.and.yarr(ipt)>ylow.and.yarr(ipt)<=yhigh) then 1259 mat(ix,iy)=mat(ix,iy)+1D0 1260 end if 1261 end do 1262 end do 1263end do 1264! !$OMP END PARALLEL DO 1265end subroutine 1266 1267 1268 1269end module 1270 1271 1272 1273 1274 1275 1276 1277!!--------- Lagrange interpolation in 1D, produce interpolated value, 1st and 2nd derivatives 1278!NOTE: 4 adjacent data points will be used to interpolate 1279!ptpos and ptval are the position and value of the input array, npt is the number of its element. ptpos must vary from small to large 1280!r is the point to be studied, the resultant val, der1, der2 are its value, 1st and 2nd derivatives 1281!itype=1: only calculate value =2: also calculate 1st-derv. =3: also calculate 2nd-derv. 1282subroutine lagintpol(ptpos,ptval,npt,r,val,der1,der2,itype) 1283implicit real*8(a-h,o-z) 1284integer npt,itype 1285real*8 ptpos(npt),ptval(npt),r,val,der1,der2 1286if (r<=ptpos(1)) then !Out of boundary 1287 der1=(ptval(2)-ptval(1))/(ptpos(2)-ptpos(1)) 1288 der2=0D0 1289 val=ptval(1)-(ptpos(1)-r)*der1 !Use linear interpolation 1290 return 1291else if (r>=ptpos(npt)) then 1292 val=0D0 !Because this function is mainly used to interpolate radial atomic density, at long distance the value must be zero 1293 der1=0D0 1294 der2=0D0 1295! val=ptval(npt) 1296! der1=(ptval(npt)-ptval(npt-1))/(ptpos(npt)-ptpos(npt-1)) 1297 return 1298end if 1299do i=1,npt !Determine which four data points will be used to interpolation 1300 if (ptpos(i)>r) exit 1301end do 1302! if (i==npt+1) i=npt !i==npt+1 means i exceeded the last point 1303iend=i+1 1304istart=i-2 1305if (istart<1) then 1306 istart=istart+1 1307 iend=iend+1 1308else if (iend>npt) then 1309 iend=iend-1 1310 istart=istart-1 1311end if 1312!Calculate interpolated value 1313val=0D0 1314do m=istart,iend 1315 poly=1D0 1316 do j=istart,iend 1317 if (j==m) cycle 1318 poly=poly* (r-ptpos(j))/(ptpos(m)-ptpos(j)) 1319 end do 1320 val=val+ptval(m)*poly 1321end do 1322!Calculate interpolated 1st-derv. 1323if (itype<2) return 1324der1=0D0 1325do m=istart,iend 1326 suml=0D0 1327 do l=istart,iend 1328 if (l==m) cycle 1329 poly=1D0 1330 do j=istart,iend 1331 if (j==m.or.j==l) cycle 1332 poly=poly* (r-ptpos(j))/(ptpos(m)-ptpos(j)) 1333 end do 1334 suml=suml+poly/(ptpos(m)-ptpos(l)) 1335 end do 1336 der1=der1+ptval(m)*suml 1337end do 1338!Calculate interpolated 2nd-derv. 1339if (itype<3) return 1340der2=0D0 1341do m=istart,iend 1342 suml=0D0 1343 do l=istart,iend 1344 if (l==m) cycle 1345 sumn=0D0 1346 do n=istart,iend 1347 if (n==l.or.n==m) cycle 1348 poly=1D0 1349 do j=istart,iend 1350 if (j==l.or.j==m.or.j==n) cycle 1351 poly=poly* (r-ptpos(j))/(ptpos(m)-ptpos(j)) 1352 end do 1353 sumn=sumn+poly/(ptpos(m)-ptpos(n)) 1354 end do 1355 suml=suml+sumn/(ptpos(m)-ptpos(l)) 1356 end do 1357 der2=der2+ptval(m)*suml 1358end do 1359end subroutine 1360