program cubcut implicit none integer*4 i,nd(3),nc,ix(3),N,nats,a,b,c,ia(2),ic real*8 r0(3),dd,v(3),dr(3),s,x(3),bohr,xx(3),der,s0, 1off character*80 fn real*8,allocatable::rh(:),r(:) logical latom c bohr=0.529177d0 write(6,6000) 6000 format(/,' cut through Gaussian .cub file',/) if(iargc().ne.1)then write(6,*)' Usage: cubcut ' stop endif call getarg(1,fn) open(9,file=fn) read(9,*) read(9,*) read(9,*)nats,(r0(i),i=1,3) read(9,*)nd(1),dr(1) read(9,*)nd(2),dr(2),dr(2) read(9,*)nd(3),dr(3),dr(3),dr(3) nc=nd(1)*nd(2)*nd(3) allocate(rh(nc),r(3*nats)) do 2 i=1,nats 2 read(9,*)(r(a+3*(i-1)),a=1,2),(r(a+3*(i-1)),a=1,3) do 3 a=1,nd(1) do 3 b=1,nd(2) 3 read(9,4000)(rh(c+(b-1)*nd(3)+(a-1)*nd(3)*nd(2)),c=1,nd(3)) 4000 format(6E13.5) close(9) do 1 a=1,3 dr(a)=dr(a)*bohr r0(a)=r0(a)*bohr do 1 i=1,nats 1 r(a+3*(i-1))=r(a+3*(i-1))*bohr write(6,*)'cube read in' write(6,600)r0(1),r0(1)+dble(nd(1))*dr(1), 1r0(2),r0(2)+dble(nd(2))*dr(2), 1r0(3),r0(3)+dble(nd(3))*dr(3),nats,nd(1),nd(2),nd(3), 1dr(1),dr(2),dr(3) 600 format(' xmin ',f12.3,' xmax ',f12.4,/, 1 ' ymin ',f12.3,' ymax ',f12.4,/, 1 ' zmin ',f12.3,' zmax ',f12.4,/, 1 i5,' atoms',/, 1 i5,' x ',i4,' x ',i4,' points ',/, 1 ' dx: ',f8.3,' dy: ',f8.3,' dz: ',f8.3,/) do 7 a=1,3 v(a)=1.0d0 7 x(a)=0.0d0 dd=0.05d0 N=100 latom=.false. off=0.0d0 ia(1)=1 ia(2)=2 44 if(latom)then call di(N,r,ia,off,v,dd,x) write(6,602)ia(1),ia(2),off,dd,N 602 format(' (1) Atoms (i1 i2) :',2i3,/, 1 ' (2) Offset :',F10.3,' A',/, 1 ' Increment :',F10.3,' A',/, 1 ' (4) Number of points :',i10,/, 1 ' (5) Definition through coordinates ',/, 1 ' (6) Continue',/, 1 ' ---------------------------------',/, 1 ' Input your choice: ',$) read(5,*)ic goto(81,82,83,84,85,86)ic goto 44 81 write(6,810) 810 format(' Atom1 atom2: ',$) read(5,*)ia(1),ia(2) goto 44 82 write(6,820) 820 format(' Offset (in A): ',$) read(5,*)off goto 44 83 goto 44 84 write(6,840) 840 format(' Number of points: ',$) read(5,*)N goto 44 85 latom=.false. goto 44 86 call di(N,r,ia,off,v,dd,x) else write(6,601)(x(a),a=1,3),(v(a),a=1,3),dd,N 601 format(' (1) Origin (ox oy oz ):',3F10.3,' A',/, 1 ' (2) Direction vector (vx vy vz):',3F10.3,/, 1 ' (3) Increment :',F10.3,' A',/, 1 ' (4) Number of points :',i10,/, 1 ' (5) Definition through atoms ',/, 1 ' (6) Continue',/, 1 ' ---------------------------------',/, 1 ' Input your choice: ',$) read(5,*)ic goto(71,72,73,74,75,76)ic goto 44 71 write(6,710) 710 format(' Initial point x y z (A): ',$) read(5,*)x(1),x(2),x(3) goto 44 72 write(6,720) 720 format(' Direction vector vx vy vz: ',$) read(5,*)v(1),v(2),v(3) goto 44 73 write(6,730) 730 format(' Increment (A): ',$) read(5,*)dd goto 44 74 write(6,740) 740 format(' Number of points: ',$) read(5,*)N goto 44 75 latom=.true. goto 44 76 call nv(v) endif open(9,file='X.PRN') do 6 i=1,N do 61 a=1,3 xx(a)=x(a)+dd*v(a)*dble(i-1) 61 ix(a)=INT((xx(a)-r0(a))/dr(a))+1 s=0.0d0 if(ix(1).ge.1.and.ix(1).le.nd(1).and. 1 ix(2).ge.1.and.ix(2).le.nd(2).and. 1 ix(3).ge.1.and.ix(3).le.nd(3))then s0=rh(ix(3)+(ix(2)-1)*nd(3)+(ix(1)-1)*nd(3)*nd(2)) s=s0 do 62 a=1,3 if(ix(a)+1.le.nd(a))then ix(a)=ix(a)+1 der=(rh(ix(3)+(ix(2)-1)*nd(3)+(ix(1)-1)*nd(3)*nd(2))-s0)/dr(a) ix(a)=ix(a)-1 s=s+der*(xx(a)-r0(a)-dr(a)*dble(ix(a)-1)) endif 62 continue else s=0.0d0 endif 6 write(9,900)dd*dble(i-1),s 900 format(f8.3,g15.5) close(9) write(6,*)'X.PRN written' goto 44 end subroutine nv(v) real*8 v(*),vn vn=dsqrt(v(1)**2+v(2)**2+v(3)**2) v(1)=v(1)/vn v(2)=v(2)/vn v(3)=v(3)/vn return end subroutine di(N,r,ia,off,v,dd,x) c define v, first point, and increment real*8 r(*),v(*),x(*),off,dd,dist,y(3) integer*4 a,ia(*),N do 1 a=1,3 1 v(a)=r(a+3*(ia(2)-1))-r(a+3*(ia(1)-1)) call nv(v) do 2 a=1,3 y(a)=r(a+3*(ia(2)-1))+off*v(a) 2 x(a)=r(a+3*(ia(1)-1))-off*v(a) dist=dsqrt((x(1)-y(1))**2+(x(2)-y(2))**2+(x(3)-y(3))**2) dd=dist/dble(N-1) return end