PROGRAM SURFACES IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) parameter (nvar0=4,ngp0=14000,nat0=100,n90=9) logical lper(nvar0),lprop dimension ngr(nvar0),ity(nat0),l1(ngp0),nh(nvar0),nbas(nvar0), 1l2(ngp0),l3(ngp0),nbase(nvar0),peri(nvar0), 2ngrh(nvar0),ice(nvar0),ic(nvar0),ich(nvar0),x0(nvar0),am0(nvar0), 3ip(nvar0),w(nvar0),l1n(ngp0),l2n(ngp0),l3n(ngp0),list(n90) real*8 min(nvar0),max(nvar0),var(ngp0,nvar0),e(ngp0), 1r(3*nat0,ngp0),a2(2,2),a2i(2,2),e2(2,4),prop(ngp0), 2x9(n90),y9(n90),dl9(n90),e12(6,12),a6(6,6),a6i(6,6), 3b6(6),v6(6),htocm real*4 amas(89) data amas/1.007825,4.003, 2 6.941, 9.012, 10.810,12.0 ,14.003074,15.994915,18.9984,20.179, 3 22.990,24.305, 26.981,28.086,30.9737,31.972,34.968852,39.948, 4 39.098,40.080,44.956,47.900,50.941,51.996,54.938,55.847, 4 58.933,58.700,63.546,65.380, 4 69.720,72.590,74.922,78.960,79.904,83.800, 5 85.468,87.620,88.906,91.220,92.906,95.940,98.906,101.070, 5 102.906,106.400,107.868,112.410, 5 114.82,118.69,121.75,127.600,126.905,131.300, 6 132.905,137.330,138.906, 6 140.120,140.908,144.240,145.000,150.400, 6 151.960,157.250,158.925,162.500,164.930,167.260,168.934, 6 173.040,174.970, 6 178.490,180.948,183.850,186.207,190.200,192.220,195.090, 6 196.967,207.590,204.370,207.200,208.980,210.000,210.001, 6 222.02, 7 223.000,226.025,227.028/ character*5 s5 c write(6,6000) 6000 format(' Molecular PESs',/,/, 1' Input : SURF.OPT ... options',/, 1' EN.LST ... list of variables and energies',/, 1' PROPERTY.LST ... list of a property (Optional)',/, 1' FILEC.X ... list of geometries',/,/, 1' Output: ENN.LST ... resampled list of geometries and G',/) C call getpar(nvar,nvar0,lper,ngr,min,max,nbase,nbas,nh,ngrh, 1ice,ic,ich,wcm,x0,am0,ipl,ip,ikin,kindpw,kinddia,T,w,99999,9999, 2peri,n9) htocm=219474.63d0 ie=1 open(8,file='EN.LST') 1 read(8,*,end=88,err=88)(var(ie,ii),ii=1,nvar),e(ie) ie=ie+1 if(ie.gt.ngp0)call report('Too many grid points') goto 1 88 ie=ie-1 close(8) write(6,*)ie,' energies calculated' inquire(file='PROPERTY.LST',exist=lprop) if(lprop)then open(8,file='PROPERTY.LST') do 7 ii=1,ie 7 read(8,*)prop(ii) close(8) write(6,*)'PROPERTY.LST read' endif open(8,FILE='FILEC.X') do 5 ii=1,ie read(8,*) read(8,*)nat if(nat.gt.nat0)call report('Too many atoms') do 5 ia=1,nat read(8,*)ity(ia),(r(ix+3*(ia-1),ii),ix=1,3) c dsqrt(M)*X in atomic units: do 5 ix=1,3 amu=1822.88d0 ff=dsqrt(amu*dble(amas(ity(ia))))/0.529177d0 5 r(ix+3*(ia-1),ii)=r(ix+3*(ia-1),ii)*ff close(8) write(6,*)ie,' geometries read' c Resample the grid: open(8,file='ENN.LST') write(s5,990)int(T) 990 format(i5) do 66 is=1,5 66 if(s5(is:is).ne.' ')goto 661 661 open(81,file=s5(is:5)//'EKT.PRN') if(lprop)open(86,file='PROPERTYN.LST') if(nvar.eq.1)then dx=(max(1)-min(1))/dble(ngr(1)) x=min(1)-dx qe=0.0d0 do 2 ix=0,ngr(1) x=x+dx call findx(i1,i2,i3,var,x,lper(1),ie,ngp0,nvar0,peri(1)) q1=var(i1,1) q2=var(i2,1) q3=var(i3,1) call findx(i1n,i2n,i3n,var,x,.false.,ie,ngp0,nvar0,peri(1)) q1n=var(i1n,1) q2n=var(i2n,1) q3n=var(i3n,1) c write(6,8009)x,q1,q2,q3,q1n,q2n,q3n c8009 format(10f8.2) if(lper(1))then tp=peri(1) tph=peri(1)/2.0d0 if(q1-x.gt.tph)q1=q1-tp if(x-q1.gt.tph)q1=q1+tp if(q2-x.gt.tph)q2=q2-tp if(x-q2.gt.tph)q2=q2+tp if(q3-x.gt.tph)q3=q3-tp if(x-q3.gt.tph)q3=q3+tp endif x1=e(i1) x2=e(i2) x3=e(i3) call abc(a,b,c,q1,q2,q3,x1,x2,x3) en=x*(a*x+b)+c if(lprop)then x1=prop(i1) x2=prop(i2) x3=prop(i3) call abc(a,b,c,q1,q2,q3,x1,x2,x3) propn=x*(a*x+b)+c write(86,808)ix,x,propn endif c s11=0.0d0 do 6 i=1,nat do 6 iax=1,3 iind=iax+3*(i-1) x1=r(iind,i1n) x2=r(iind,i2n) x3=r(iind,i3n) call abc(a,b,c,q1n,q2n,q3n,x1,x2,x3) dxid1=2.0d0*x*a+b 6 s11=s11+dxid1**2 qe=qe+exp(-en*627.5d0/1.878d-3/T) 8081 format(2f12.6) 2 write(8,808)ix,x,en,1.0d0/s11 c 1/s11 ... inverted metric tensor 808 format(i5,2f12.6,2g12.5) rewind 8 do 202 ix=0,ngr(1) read(8,*)idum,x,en 202 write(81,8081)x,exp(-en*627.5d0/1.878d-3/T)/dx/qe endif c 2D============================================================== if(nvar.eq.2)then if(ice(1).eq.-2)then c c quadratic fit from 9 points: c --------------------------------------------------------------- dx=(max(1)-min(1))/dble(ngr(1)) dy=(max(2)-min(2))/dble(ngr(2)) x=min(1)-dx do 399 ix=0,ngr(1) x=x+dx y=min(2)-dy do 399 iy=0,ngr(2) write(6,6001)ix,iy,x,y 6001 format(2i4,5f15.4) y=y+dy c find 9 closest points to (x,y): call find9(x,y,lper(1),lper(2),ngp0,nvar0,ie, 1 list,9,x9,y9,dl9,peri(1),peri(2),var) c c fit quadratic spline: do i=1,6 do j=1,6 a6(i,j)=0.0d0 do ii=1,n9 f1=1.0d0 if(i.eq.2)f1=x9(ii)-x if(i.eq.3)f1=y9(ii)-y if(i.eq.4)f1=(x9(ii)-x)**2/2.0d0 if(i.eq.5)f1=(y9(ii)-y)**2/2.0d0 if(i.eq.6)f1=(y9(ii)-y)*(x9(ii)-x) f2=1.0d0 if(j.eq.2)f2=x9(ii)-x if(j.eq.3)f2=y9(ii)-y if(j.eq.4)f2=(x9(ii)-x)**2/2.0d0 if(j.eq.5)f2=(y9(ii)-y)**2/2.0d0 if(j.eq.6)f2=(y9(ii)-y)*(x9(ii)-x) a6(i,j)=a6(i,j)+f1*f2 enddo enddo enddo do i=1,6 b6(i)=0.0d0 enddo do ii=1,n9 b6(1)=b6(1)+e(list(ii)) b6(2)=b6(2)+e(list(ii))*(x9(ii)-x) b6(3)=b6(3)+e(list(ii))*(y9(ii)-y) b6(4)=b6(4)+e(list(ii))*(x9(ii)-x)*(x9(ii)-x)/2.0d0 b6(5)=b6(5)+e(list(ii))*(y9(ii)-y)*(y9(ii)-y)/2.0d0 b6(6)=b6(6)+e(list(ii))*(x9(ii)-x)*(y9(ii)-y) enddo TOL=1.0d-9 call INV(6,a6,a6i,6,TOL,e12,IERR) do i=1,6 v6(i)=0.0d0 do ii=1,6 v6(i)=v6(i)+a6i(i,ii)*b6(ii) enddo enddo c v6 components:E,dE/dx,dE/dy,d^2E/dxx,dE^2/dyy,d^2E/dxy en=v6(1) if(lprop)then do i=1,6 b6(i)=0.0d0 enddo do ii=1,n9 b6(1)=b6(1)+prop(list(ii)) b6(2)=b6(2)+prop(list(ii))*(x9(ii)-x) b6(3)=b6(3)+prop(list(ii))*(y9(ii)-y) b6(4)=b6(4)+prop(list(ii))*(x9(ii)-x)*(x9(ii)-x)/2.0d0 b6(5)=b6(5)+prop(list(ii))*(y9(ii)-y)*(y9(ii)-y)/2.0d0 b6(6)=b6(6)+prop(list(ii))*(x9(ii)-x)*(y9(ii)-y) enddo propn=0.0d0 do ii=1,6 propn=propn+a6i(i,ii)*b6(ii) enddo write(86,819)ix,iy,x,y,propn endif s11=0.0d0 s22=0.0d0 s12=0.0d0 do 499 i=1,nat do 499 iax=1,3 iind=iax+3*(i-1) do ii=1,6 b6(ii)=0.0d0 enddo do ii=1,n9 b6(1)=b6(1)+r(iind,list(ii)) b6(2)=b6(2)+r(iind,list(ii))*(x9(ii)-x) b6(3)=b6(3)+r(iind,list(ii))*(y9(ii)-y) b6(4)=b6(4)+r(iind,list(ii))*(x9(ii)-x)*(x9(ii)-x)/2.0d0 b6(5)=b6(5)+r(iind,list(ii))*(y9(ii)-y)*(y9(ii)-y)/2.0d0 b6(6)=b6(6)+r(iind,list(ii))*(x9(ii)-x)*(y9(ii)-y) enddo do jj=1,6 v6(jj)=0.0d0 do ii=1,6 v6(jj)=v6(jj)+a6i(jj,ii)*b6(ii) enddo enddo d1=v6(2) d2=v6(3) s11=s11+d1*d1 s12=s12+d1*d2 499 s22=s22+d2*d2 a2(1,1)=s11 a2(1,2)=s12 a2(2,1)=s12 a2(2,2)=s22 TOL=1.0d-9 call INV(2,a2,a2i,2,TOL,e2,IERR) det=a2i(1,1)*a2i(2,2)-a2i(1,2)*a2i(2,1) if(IERR.ne.0)call report('inversion unsuccessfull') write(81,809)x,y,exp(-en*627.5d0/1.878d-3/T) 399 write(8,819)ix,iy,x,y,en,det,a2i(1,1),a2i(1,2),a2i(2,2) c --------------------------------------------------------------- else c --------------------------------------------------------------- dx=(max(1)-min(1))/dble(ngr(1)) dy=(max(2)-min(2))/dble(ngr(2)) x=min(1)-dx do 3 ix=0,ngr(1) x=x+dx c c if symmetry, use it for energy call findlist(l1,l2,l3,var,x,lper(1),ie,ngp0,nvar0,peri(1)) q1=var(l1(2),1) q2=var(l2(2),1) q3=var(l3(2),1) if(lper(1))then tp=peri(1) tph=peri(1)/2.0d0 if(q1-x.gt.tph)q1=q1-tp if(x-q1.gt.tph)q1=q1+tp if(q2-x.gt.tph)q2=q2-tp if(x-q2.gt.tph)q2=q2+tp if(q3-x.gt.tph)q3=q3-tp if(x-q3.gt.tph)q3=q3+tp endif c c non-symmetrical points: call findlist(l1n,l2n,l3n,var,x,.false.,ie,ngp0,nvar0,peri(1)) q1n=var(l1n(2),1) q2n=var(l2n(2),1) q3n=var(l3n(2),1) c c l1 l2 l3 c c | c | c | p13 c p11 p12 | c p22 | c y-------------p21-----------+---------p23---- f1 f2 f3 values c | from Ps for q1 q2 q3 c p32 | c p31 | p33 c | c | c | c | c q1 q2 | q3 closest calculated points c | c | c x x value c y=min(2)-dy do 3 iy=0,ngr(2) y=y+dy c find 3 best ys in l1, determine value f1 call findy(y,l1,i11,i21,i31,var,lper(2),ngp0,nvar0,peri(2)) p11=var(l1(i11),2) p21=var(l1(i21),2) p31=var(l1(i31),2) c find 3 best ys in l2, determine value f2 call findy(y,l2,i12,i22,i32,var,lper(2),ngp0,nvar0,peri(2)) p12=var(l2(i12),2) p22=var(l2(i22),2) p32=var(l2(i32),2) c find 3 best ys in l3, determine value f3 call findy(y,l3,i13,i23,i33,var,lper(2),ngp0,nvar0,peri(2)) p13=var(l3(i13),2) p23=var(l3(i23),2) p33=var(l3(i33),2) if(lper(2))then tp=peri(2) tph=peri(2)/2.0d0 if(p11-y.gt.tph)p11=p11-tp if(y-p11.gt.tph)p11=p11+tp if(p21-y.gt.tph)p21=p21-tp if(y-p21.gt.tph)p21=p21+tp if(p31-y.gt.tph)p31=p31-tp if(y-p31.gt.tph)p31=p31+tp if(p12-y.gt.tph)p12=p12-tp if(y-p12.gt.tph)p12=p12+tp if(p22-y.gt.tph)p22=p22-tp if(y-p22.gt.tph)p22=p22+tp if(p32-y.gt.tph)p32=p32-tp if(y-p32.gt.tph)p32=p32+tp if(p13-y.gt.tph)p13=p13-tp if(y-p13.gt.tph)p13=p13+tp if(p23-y.gt.tph)p23=p23-tp if(y-p23.gt.tph)p23=p23+tp if(p33-y.gt.tph)p33=p33-tp if(y-p33.gt.tph)p33=p33+tp endif c c non-symmetrical call findy(y,l1n,i11n,i21n,i31n,var,.false.,ngp0,nvar0,peri(2)) p11n=var(l1n(i11n),2) p21n=var(l1n(i21n),2) p31n=var(l1n(i31n),2) call findy(y,l2n,i12n,i22n,i32n,var,.false.,ngp0,nvar0,peri(2)) p12n=var(l2n(i12n),2) p22n=var(l2n(i22n),2) p32n=var(l2n(i32n),2) call findy(y,l3n,i13n,i23n,i33n,var,.false.,ngp0,nvar0,peri(2)) p13n=var(l3n(i13n),2) p23n=var(l3n(i23n),2) p33n=var(l3n(i33n),2) g1=e(l1(i11)) g2=e(l1(i21)) g3=e(l1(i31)) call abc(a,b,c,p11,p21,p31,g1,g2,g3) f1=y*(a*y+b)+c g1=e(l2(i12)) g2=e(l2(i22)) g3=e(l2(i32)) call abc(a,b,c,p12,p22,p32,g1,g2,g3) f2=y*(a*y+b)+c g1=e(l3(i13)) g2=e(l3(i23)) g3=e(l3(i33)) call abc(a,b,c,p13,p23,p33,g1,g2,g3) f3=y*(a*y+b)+c call abc(a,b,c,q1,q2,q3,f1,f2,f3) en=x*(a*x+b)+c if(lprop)then g1=prop(l1(i11)) g2=prop(l1(i21)) g3=prop(l1(i31)) call abc(a,b,c,p11,p21,p31,g1,g2,g3) f1=y*(a*y+b)+c g1=prop(l2(i12)) g2=prop(l2(i22)) g3=prop(l2(i32)) call abc(a,b,c,p12,p22,p32,g1,g2,g3) f2=y*(a*y+b)+c g1=prop(l3(i13)) g2=prop(l3(i23)) g3=prop(l3(i33)) call abc(a,b,c,p13,p23,p33,g1,g2,g3) f3=y*(a*y+b)+c call abc(a,b,c,q1,q2,q3,f1,f2,f3) propn=x*(a*x+b)+c write(86,819)ix,iy,x,y,propn endif s11=0.0d0 s22=0.0d0 s12=0.0d0 do 4 i=1,nat do 4 iax=1,3 iind=iax+3*(i-1) g1=r(iind,l1n(i11n)) g2=r(iind,l1n(i21n)) g3=r(iind,l1n(i31n)) call abc(a,b,c,p11n,p21n,p31n,g1,g2,g3) f1=2.0d0*y*a+b r1=y*(a*y+b)+c g1=r(iind,l2n(i12n)) g2=r(iind,l2n(i22n)) g3=r(iind,l2n(i32n)) call abc(a,b,c,p12n,p22n,p32n,g1,g2,g3) f2=2.0d0*y*a+b r2=y*(a*y+b)+c g1=r(iind,l3n(i13n)) g2=r(iind,l3n(i23n)) g3=r(iind,l3n(i33n)) call abc(a,b,c,p13n,p23n,p33n,g1,g2,g3) f3=2.0d0*y*a+b r3=y*(a*y+b)+c call abc(a,b,c,q1n,q2n,q3n,f1,f2,f3) d2=x*(a*x+b)+c call abc(a,b,c,q1n,q2n,q3n,r1,r2,r3) d1=2.0d0*x*a+b s11=s11+d1*d1 s12=s12+d1*d2 4 s22=s22+d2*d2 a2(1,1)=s11 a2(1,2)=s12 a2(2,1)=s12 a2(2,2)=s22 TOL=1.0d-9 call INV(2,a2,a2i,2,TOL,e2,IERR) det=a2i(1,1)*a2i(2,2)-a2i(1,2)*a2i(2,1) if(IERR.ne.0)call report('inversion unsuccessfull') write(81,809)x,y,exp(-en*627.5d0/1.878d-3/T) 809 format(2f12.6,g13.5) 3 write(8,819)ix,iy,x,y,en,det,a2i(1,1),a2i(1,2),a2i(2,2) 819 format(2i5,3f12.6,4g13.5) endif c --------------------------------------------------------------- endif c ============================================================== close(8) close(86) close(81) stop end c ------------------------------------------------------------ subroutine report(s) character*(*) s write(6,*)s stop end c ------------------------------------------------------------ subroutine findx(i1,i2,i3,var,x,lp1,ie,ngp0,nvar0,tp) c finds 3 y-arrays in x1 x2 x3 closest to x IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) logical lp1,nonequal real*8 var(ngp0,nvar0) to=1.0d-4 c c initialize x1 as var(1 i1=1 d1=ddi(lp1,var(i1,1),x,tp) c c get some x2 different from x1 do 10 i2=2,ie 10 if(nonequal(lp1,var(i2,1),var(i1,1),to,tp))goto 11 call report('x2 not found in findx') 11 d2=ddi(lp1,var(i2,1),x,tp) c c get x3 different from x1 and x2 do 12 i3=3,ie 12 if(nonequal(lp1,var(i3,1),var(i1,1),to,tp).and. 1nonequal(lp1,var(i3,1),var(i2,1),to,tp))goto 13 call report('x3 not found in findx') 13 d3=ddi(lp1,var(i3,1),x,tp) call order(d1,d2,d3,i1,i2,i3) c look if some points are closer: do 1 ii=4,ie dist=ddi(lp1,var(ii,1),x,tp) if(dist.lt.d1-to)then d3=d2 i3=i2 d2=d1 i2=i1 d1=dist i1=ii endif if(dist.lt.d2-to.and.nonequal(lp1,var(ii,1),var(i1,1),to,tp))then d3=d2 i3=i2 d2=dist i2=ii endif if(dist.lt.d3-to.and.nonequal(lp1,var(ii,1),var(i1,1),to,tp).and. 1nonequal(lp1,var(ii,1),var(i2,1),to,tp))then d3=dist i3=ii endif 1 continue return end c ------------------------------------------------------------ subroutine findlist(l1,l2,l3,var,x,lp1,ie,ngp0,nvar0,tp) c finds 3 y-arrays in x1 x2 x3 closest to x IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) logical lp1,nonequal real*8 var(ngp0,nvar0) dimension l1(*),l2(*),l3(*) to=1.0d-4 call findx(i1,i2,i3,var,x,lp1,ie,ngp0,nvar0,tp) c record list of ys for xs same as x1,x2,x3 nl1=1 nl2=1 nl3=1 do 2 ii=1,ie if(.not.nonequal(lp1,var(ii,1),var(i1,1),to,tp))then nl1=nl1+1 l1(nl1)=ii endif if(.not.nonequal(lp1,var(ii,1),var(i2,1),to,tp))then nl2=nl2+1 l2(nl2)=ii endif if(.not.nonequal(lp1,var(ii,1),var(i3,1),to,tp))then nl3=nl3+1 l3(nl3)=ii endif 2 continue l1(1)=nl1 l2(1)=nl2 l3(1)=nl3 return end c ------------------------------------------------------------ subroutine findy(y,l1,i1,i2,i3,var,lp,ngp0,nvar0,tp) c in l1 find y1 y2 y3 closest to y IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) logical lp,nonequal real*8 var(ngp0,nvar0) dimension l1(*) to=1.0d-4 c initalize y1 i1=2 d1=ddi(lp,var(l1(i1),2),y,tp) c find y2 different from y1 do 10 i2=3,l1(1) 10 if(nonequal(lp,var(l1(i2),2),var(l1(i1),2),to,tp))goto 11 write(6,*)' y =',y,' tol =',to write(6,*)' tp =',tp do 101 i2=3,l1(1) 101 write(6,*)i2,lp,var(l1(i2),2),var(l1(i1),2) call report('y2 not found in findy') 11 d2=ddi(lp,var(l1(i2),2),y,tp) c find y3 different from y1 and from y2 do 12 i3=4,l1(1) 12 if(nonequal(lp,var(l1(i3),2),var(l1(i1),2),to,tp).and. 1nonequal(lp,var(l1(i3),2),var(l1(i2),2),to,tp))goto 13 call report('y3 not found in findy') 13 d3=ddi(lp,var(l1(i3),2),y,tp) call order(d1,d2,d3,i1,i2,i3) c look for closer points: do 1 ii=5,l1(1) dist=ddi(lp,var(l1(ii),2),y,tp) if(dist.lt.d1-to)then d3=d2 i3=i2 d2=d1 i2=i1 d1=dist i1=ii endif if(dist.lt.d2-to.and.nonequal(lp,var(l1(ii),2),var(l1(i1),2),to,tp 1))then d3=d2 i3=i2 d2=dist i2=ii endif if(dist.lt.d3-to.and.nonequal(lp,var(l1(ii),2),var(l1(i1),2),to,tp 1).and.nonequal(lp,var(l1(ii),2),var(l1(i2),2),to,tp))then d3=dist i3=ii endif 1 continue return end c --------------------------------------------------------------- function ddi(l,x,y,tp) c distance between two coordinates, including periodicity real*8 x,y,tp,d,ddi logical l d=abs(x-y) if(l)then if(d.gt.abs(x-y-tp))d=abs(x-y-tp) if(d.gt.abs(x-y+tp))d=abs(x-y+tp) endif ddi=d return end c ------------------------------------------------------------ function nonequal(l,x,y,to,tp) logical nonequal,l real*8 x,y,to,tp if((l.and.abs(x-y).gt.to.and.abs(x-y+tp).gt.to.and. 1abs(x-y-tp).gt.to).or.(.not.l.and.abs(x-y).gt.to))then nonequal=.true. else nonequal=.false. endif return end c ------------------------------------------------------------ subroutine abc(a,b,c,q1,q2,q3,x1,x2,x3) real*8 a,b,c,q1,q2,q3,x1,x2,x3,down down=-(q1-q2)*(q2-q3)*(q3-q1) a=(x1*(q2-q3)+x2*(q3-q1)+x3*(q1-q2))/down b=-(x1*(q2**2-q3**2)+x2*(q3**2-q1**2)+x3*(q1**2-q2**2))/down c=(x1*(q2-q3)*q2*q3+x2*(q3-q1)*q1*q3+x3*(q1-q2)*q1*q2)/down return end c ------------------------------------------------------------ subroutine order(d1,d2,d3,i1,i2,i3) real*8 d1,d2,d3,t integer*4 i1,i2,i3 if(d1.gt.d2)then t=d1 d1=d2 d2=t it=i1 i1=i2 i2=it endif if(d1.gt.d3)then t=d1 d1=d3 d3=t it=i1 i1=i3 i3=it endif if(d2.gt.d3)then t=d2 d2=d3 d3=t it=i2 i2=i3 i3=it endif return end c =========================================== SUBROUTINE INV(nr,a,ai,n,TOL,e,IERR) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) dimension ai(nr,nr),a(nr,nr),e(nr,2*nr) C 10000 IERR=0 DO 1 ii=1,n DO 1 jj=1,n e(ii,jj)=a(ii,jj) 1 e(ii,jj+n)=0.0D0 do 13 ii=1,n 13 e(ii,ii+n)=1.0D0 c DO 2 ii=1,n-1 if (ABS(e(ii,ii)).LE.TOL) then DO 3 io=ii+1,n 3 if (ABS(e(io,ii)).GT.TOL) goto 11 IERR=1 write(6,*)ii write(6,*)'tol = ',tol tol=tol*0.50d0 if(tol.gt.1.0d-20)goto 10000 RETURN c 11 CONTINUE DO 4 kk=1,2*n w=e(ii,kk) e(ii,kk)=e(io,kk) 4 e(io,kk)=w ENDIF eii=e(ii,ii) DO 5 jj=ii+1,n e1=e(jj,ii)/eii DO 6 kk=ii+1, 2*n 6 e(jj,kk)=e(jj,kk)-e(ii,kk)*e1 5 e(jj,ii)=0.0D0 2 CONTINUE c DO 7 i2=n,2,-1 eii=e(i2,i2) DO 7 j2=i2-1,1,-1 e1=e(j2,i2)/eii DO 9 kk=1, n 9 e(j2,kk+n)=e(j2,kk+n)-e(i2,kk+n)*e1 7 e(j2,i2)=0.0d0 c DO 10 ii=1,n ei=1.0d0/e(ii,ii) DO 12 jj=1,n 12 ai(ii,jj)=e(ii,jj+n)*ei 10 CONTINUE c RETURN END subroutine getpar(nvar,nvar0,lper,ngr,min,max,nbase,nbas,nh,ngrh, 1ice,ic,ich,wcm,x0,am0,ipl,ip,ikin,kindpw,kinddia,T,w,nbas0,nh0, 2peri,n9) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) character*4 key c icon=0 ... reads par for surfaces c icon=1 ... reads par for peswf dimension ip(*),ngr(*),min(*),max(*),ngrh(*),nbase(*), 1nbas(*),nh(*),ice(*),ic(*),ich(*),x0(*),am0(*),w(*),lper(*), 2peri(*) logical lper real*8 min,max n9=9 ipl=0 ikin=1 kindpw=1 kinddia=1 T=298.15d0 il=0 open(8,file='SURF.OPT') 1 read(8,800,end=88,err=88)key 800 format(a4) il=il+1 c c Obligatory Parameters: if(key.eq.'NVAR'.and.il.eq.1)then read(8,*)nvar if(nvar.gt.nvar0)call report('Too many variables') do 2 ii=1,nvar write(6,6005)ii 6005 format(' Variable',i2,':') peri(ii)=360.0d0 read(8,*)lper(ii) read(8,*)ngr(ii) read(8,*)min(ii),max(ii) read(8,*)ngrh(ii) read(8,*)nbase(ii),nbas(ii),nh(ii) read(8,*)ice(ii),ic(ii),ich(ii) read(8,*)wcm,x0(ii),am0(ii) if(ice(ii).eq.0)then write(6,6002)nbase(ii), ' plane waves','E' 6002 format(i4,a12,' for ',a1) else if(ice(ii).eq.1)then write(6,6002)nbase(ii),' harmonics','E' else if(ice(ii).eq.-1)then write(6,*)'Quadratic fit for E' else if(ice(ii).eq.-2)then write(6,*)'Quadratic fit for E from closest 9 points in 2D' else call report('unknown fit required') endif endif endif endif if(ic(ii).eq.0)then write(6,6002)nbas(ii), ' plane waves','G' else if(ic(ii).eq.1)then write(6,6002)nbas(ii),' harmonics','G' else write(6,*)'Quadratic fit for G' endif endif if(ich(ii).eq.0)then write(6,6002)nh(ii), ' plane waves','H' else if(ich(ii).eq.1)then write(6,6002)nh(ii),' harmonics','H' else call report('Unknown basis option for H') endif endif if(nbase(ii).gt.nbas0)call report('too many E basis functions') if(nbas(ii). gt.nbas0)call report('too many G basis functions') if(nh(ii). gt.nh0 )call report('too many H basis functions') if(ic(ii).eq.1.or.ich(ii).eq.1)write(6,5006)wcm,x0(ii),am0(ii) 5006 format(' HO basis',f12.2,' cm-1, x0 = ',f12.6,' mass = ',f10.2) 2 w(ii)=wcm/htocm else if(il.eq.1)call report('NVAR must be defined at the first line.') endif c c Optional Parameters if(key.eq.'STAT')read(8,*)ipl,(ip(ii),ii=1,ipl) if(key.eq.'KINE')read(8,*)ikin if(key.eq.'PLAN')read(8,*)kindpw if(key.eq.'DIAG')read(8,*)kinddia if(key.eq.'TEMP')read(8,*)T if(key.eq.'PERI')read(8,*) (peri(i),i=1,NVAR) goto 1 88 close(8) return end c ------------------------------------------------------------ subroutine find9(x,y,lp1,lp2,ngp0,nvar0,ie,list,n9,x9,y9,dl9, 1tp1,tp2,var) c finds 9 points closest to (x,y) implicit none integer*4 ngp0,nvar0,i,ie,list(*),n9 logical lp1,lp2 real*8 var(ngp0,nvar0),dist,x9(*),y9(*),dl9(*),tp1,tp2,xe,ye, 1x,y c c take first 9 points do 1 i=1,n9 list(i)=i x9(i)=var(i,1) y9(i)=var(i,2) dl9(i)=(x-var(i,1))**2+(y-var(i,2))**2 if(lp1)then if((x+tp1-var(i,1))**2+(y-var(i,2))**2.lt.dl9(i))then x9(i)=var(i,1)-tp1 dl9(i)=(x+tp1-var(i,1))**2+(y-var(i,2))**2 endif if((x-tp1-var(i,1))**2+(y-var(i,2))**2.lt.dl9(i))then x9(i)=var(i,1)+tp1 dl9(i)=(x-tp1-var(i,1))**2+(y-var(i,2))**2 endif endif if(lp2)then if((x-var(i,1))**2+(y-tp2-var(i,2))**2.lt.dl9(i))then y9(i)=var(i,2)+tp2 dl9(i)=(x-var(i,1))**2+(y-tp2-var(i,2))**2 endif if((x-var(i,1))**2+(y+tp2-var(i,2))**2.lt.dl9(i))then y9(i)=var(i,2)-tp2 dl9(i)=(x-var(i,1))**2+(y+tp2-var(i,2))**2 endif endif 1 continue call order9(list,x9,y9,dl9,n9) c c go through rest of points, look if they are better do 2 i=n9+1,ie xe=var(i,1) ye=var(i,2) dist=(x-xe)**2+(y-ye)**2 if(lp1)then if((x+tp1-xe)**2+(y-ye)**2.lt.dist)then xe=xe-tp1 dist=(x-xe)**2+(y-ye)**2 endif if((x-tp1-xe)**2+(y-ye)**2.lt.dist)then xe=xe+tp1 dist=(x-xe)**2+(y-ye)**2 endif endif if(lp2)then if((x-xe)**2+(y-tp2-ye)**2.lt.dist)then ye=ye+tp2 dist=(x-xe)**2+(y-ye)**2 endif if((x-xe)**2+(y+tp2-ye)**2.lt.dist)then ye=ye-tp2 dist=(x-xe)**2+(y-ye)**2 endif endif if(dist.lt.dl9(n9))then dl9(n9)=dist list(n9)=i x9(n9)=xe y9(n9)=ye call order9(list,x9,y9,dl9,n9) endif 2 continue return end c ------------------------------------------------------------ subroutine order9(list,x9,y9,dl9,n) implicit none integer*4 it,list(*),i,j,n real*8 t,x9(*),y9(*),dl9(*) do 1 i=1,n do 1 j=i+1,n if(dl9(j).lt.dl9(i))then t=dl9(i) dl9(i)=dl9(j) dl9(j)=t t=x9(i) x9(i)=x9(j) x9(j)=t t=y9(i) y9(i)=y9(j) y9(j)=t it=list(i) list(i)=list(j) list(j)=it endif 1 continue return end