program ahpes c Sigmaplot AH surface plotting from AH constants implicit none integer*4 nq,n,nat,np,i,j,ii,jj real*8 qmin,qmax,dq,qi,qj,V,Vmax,Vmin,a(2,2),CM logical lex integer*4,allocatable::z(:) real*8,allocatable::smat(:,:),w(:),c(:,:,:),d(:,:,:),p(:,:,:), 1s(:,:,:) CM=219474.0d0 write(6,900) 900 format(' ',/,/, 1 ' Input:',/, 1 ' F.INP old S-matrix',/, 1 ' CQQ.SCR.TXT cubic force constants',/, 1 ' DQQ.SCR.TXT quartic force constants',/, 1 ' 5QQ.SCR.TXT 5 force constants',/, 1 ' 6QQ.SCR.TXT 6 force constants',/, 1 ' Output:',/, 1 ' PES.TXT PES profile',/) open(9,file='F.INP',status='old') read(9,*)nq,n,nat close(9) write(6,*)nq,' modes, ',nat,' atoms in F.INP' write(6,*)'Modes to scan I J:' read(5,*)i,j write(6,*)'Qmin Qmax Np:' read(5,*)qmin,qmax,NP allocate(z(nat),w(n),smat(n,n),c(n,n,n),d(n,n,n),p(n,n,n), 1s(n,n,n)) call reads(n,smat,w,'F.INP',z) call load3(c,n,w) call load4(d,n,w,4) call load4(p,n,w,5) call load4(s,n,w,6) c Vij = wi^2/2 qi^2 + wj^2/2 qj^2 + (1/6) (ciii+cjjj+3 cijj) c +(1/24) (diiii+djjjj+4diiij+6diijj (iijj ijij ijji jiij jiji jjii) c +(1/120)(piiiii+pjjjjj+diiij+6diijj (iijj ijij ijji jiij jiji jjii) c etc Vmax=0.0d0 Vmin=0.0d0 open(9,file='PES.TXT') dq=(qmax-qmin)/dble(Np) qi=qmin-dq do 1 ii=1,Np qi=qi+dq qj=qmin-dq do 1 jj=1,Np qj=qj+dq V=1.0d0/ 2.0d0*(w(i)*qi**2+w(j)*qj**2) 1+ 1.0d0/ 6.0d0*(c(i,i,i)*qi**3 +c(j,j,j)*qj**3) 1+ 3.0d0/ 6.0d0*(c(i,i,j)*qi*qi*qj+c(j,j,i)*qj*qj*qi) 1+ 1.0d0/ 24.0d0*(d(i,i,i)*qi**4 +d(j,j,j)*qj**4) 1+ 4.0d0/ 24.0d0*(d(i,i,j)*qi**3*qj+d(j,j,i)*qj**3*qi) 1+ 6.0d0/ 24.0d0*d(i,j,j)*qi**2*qj**2 1+ 1.0d0/120.0d0*(p(i,i,i)*qi**5+p(j,j,j)*qj**5) 1+ 5.0d0/120.0d0*(p(i,i,j)*qi**4*qj+p(j,j,i)*qj**4*qi) 1+10.0d0/120.0d0*(p(i,j,j)*qi**3*qj**2+p(j,i,i)*qj**3*qi**2) 1+ 1.0d0/720.0d0*(s(i,i,i)*qi**6 +s(j,j,j)*qj**6) 1+16.0d0/720.0d0*(s(i,i,j)*qi**5*qj +s(j,j,i)*qj**5*qi) 1+15.0d0/720.0d0*(s(i,j,j)*qi**4*qj**2+s(j,i,i)*qj**4*qi**2) if(V.gt.Vmax)Vmax=V if(V.lt.Vmin)Vmin=V 1 write(9,90)qi,qj,v 90 format(2f12.6,f12.2) close(9) write(6,65)Vmax,Vmin 65 format(' Vmax Vmin ',2f12.2) open(9,file='PESH.TXT') qi=qmin-dq do 2 ii=1,Np qi=qi+dq qj=qmin-dq do 2 jj=1,Np qj=qj+dq V=1.0d0/ 2.0d0*(w(i)*qi**2+w(j)*qj**2) 2 write(9,90)qi,qj,v close(9) inquire(file='AIJ.SCR.TXT',exist=lex) if(lex)then open(9,file='AIJ.SCR.TXT') read(9,*)i,j,a(1,1) read(9,*)i,j,a(1,2) read(9,*)i,j,a(2,2) close(9) open(9,file='PEAIJ.TXT') qi=qmin-dq do 3 ii=1,Np qi=qi+dq qj=qmin-dq do 3 jj=1,Np qj=qj+dq V=0.5d0*(a(1,1)*qi**2+a(2,2)*qj**2+2.0d0*a(1,2)*qi*qj) 3 write(9,90)qi,qj,v*CM close(9) endif end SUBROUTINE READS(N3,S,E,fn,z) IMPLICIT none integer*4 N3,z(*),i,ix,j real*8 S(N3,N3),E(*) character*(*)fn OPEN(4,FILE=fn,STATUS='OLD') read(4,*) do 1 i=1,N3/3 1 read(4,*)z(i) read(4,*) DO 2 I=1,N3/3 DO 2 J=1,N3 2 read(4,*)(s(3*(i-1)+ix,N3-j+1),ix=1,2), 1(s(3*(i-1)+ix,N3-j+1),ix=1,3) read(4,*) READ(4,4000)(E(N3-I+1),I=1,N3) 4000 FORMAT(6F11.6) close(4) write(6,*)fn//' read' RETURN end subroutine load3(C,N,e) implicit none real*8 C(N,N,N),AC,e(*) integer*4 I,J,K,N,nqt call tz(c,N) nqt=0 OPEN(4,file='CQQ.SCR.TXT') 8 READ(4,*,end=88,err=88)I,J,K,AC,AC C(I,J,K)=AC C(I,K,J)=AC C(J,K,I)=AC C(J,I,K)=AC C(K,J,I)=AC C(K,I,J)=AC nqt=nqt+1 goto 8 88 CLOSE(4) write(6,881)nqt 881 format(I10,' constants read from CQQ.SCR.TXT') return end subroutine load4(d,n,e,ic) implicit none integer*4 n,I,J,K,IT,ic real*8 d(n,n,n),AC,e(*) call tz(d,N) IT=0 if(ic.eq.4)open(23,file='DQQ.SCR.TXT') if(ic.eq.5)open(23,file='5QQ.SCR.TXT') if(ic.eq.6)open(23,file='6QQ.SCR.TXT') 7 read(23,*,err=7000,end=7000)I,J,K,AC,AC IT=IT+1 d(I,J,K)=AC d(I,K,J)=AC goto 7 7000 close(23) if(ic.eq.4)WRITE(6,7004)IT,'quartics' if(ic.eq.5)WRITE(6,7004)IT,'quintics' if(ic.eq.6)WRITE(6,7004)IT,'sextics ' 7004 FORMAT(I10,' ',a8,' read') return end subroutine tz(d,N) implicit none integer*4 I,J,K,N real*8 d(N,N,N) do 8 I=1,N do 8 J=1,N do 8 K=1,N 8 d(I,J,K)=0.0d0 return end