PROGRAM scaleam implicit none character*80 s80 integer*4 n,nat,iargc,i,ib,ix real*8 sf(5) real*8,allocatable::r(:),F(:,:) integer*4,allocatable::ic(:),io(:),kb(:) integer*4,allocatable::z(:) write(6,600) 600 format(' Scales C=O force field',/,/, 1 ' Input: FILE.FC force field',/ 1 ' FILE.X geometry',/ 1 ' Output: FILES.FC scaled force field',/,/, 1 ' SF order: 1 unknown ',/, 1 ' 2 -NH-CO-C ',/, 1 ' 3 NH2-CO-C ',/, 1 ' 4 =N-CO-C ',/, 1 ' 5 O-CO-C ',/,/) if(iargc().lt.1)call report('Usage: scaleam ..') do 2 i=1,5 if(iargc().ge.i)then call getarg(i,s80) read(s80,*)sf(i) else sf(i)=1.0d0 endif 2 continue open(9,file='FILE.X') read(9,*) read(9,*)nat n=3*nat allocate(r(n),z(nat)) do 1 i=1,nat 1 read(9,*)z(i),(r(ix+3*(i-1)),ix=1,3) close(9) write(6,*)nat,' atoms in FILE.X' allocate(F(N,N)) CALL READFF(N,F) allocate(ic(n),io(n),kb(n)) call defa(nat,r,z,ib,ic,io,kb,sf) call scale(ib,ic,io,r,f,sf,n,kb) CALL WRITEFF(N,F) end subroutine scale(ib,ic,io,r,f,sf,n,kb) implicit none integer*4 i,ic(*),io(*),ind(6),ix,iy,j,k,kk,ia,ib,n,kb(*) real*8 r(*),xc,yc,zc,xo,yo,zo,d,s(6,6),sp,spro,fi(6,6), 1f(n,n),sf(*),f6(6,6),fis(6,6),s6(6),m(6),fold,fnew write(6,*)'Old and new force constants:' do 9 ia=1,6 9 s6(ia)=1.0d0 m(1)=6.0d0 m(2)=6.0d0 m(3)=6.0d0 m(4)=8.0d0 m(5)=8.0d0 m(6)=8.0d0 do 1 i=1,ib s6(1)=sf(kb(i)) ind(1)=1+3*(ic(i)-1) ind(2)=2+3*(ic(i)-1) ind(3)=3+3*(ic(i)-1) ind(4)=1+3*(io(i)-1) ind(5)=2+3*(io(i)-1) ind(6)=3+3*(io(i)-1) xc=r(ind(1)) yc=r(ind(2)) zc=r(ind(3)) xo=r(ind(4)) yo=r(ind(5)) zo=r(ind(6)) d=dsqrt((xc-xo)**2+(yc-yo)**2+(zc-zo)**2) c define mini-S-matrix do 2 ix=1,6 do 2 iy=1,6 2 s(ix,iy)=0.0d0 do 3 iy=1,6 3 s(iy,iy)=1.0d0 c replace xc by the bond length s(1,1)= (xc-xo)/d s(2,1)= (yc-yo)/d s(3,1)= (zc-zo)/d s(4,1)=-(xc-xo)/d s(5,1)=-(yc-yo)/d s(6,1)=-(zc-zo)/d c orthonormalize sp=spro(1,1,s) call norm(s,1,sp) do 4 j=2,6 do 5 k=1,j-1 sp=spro(j,k,s) do 5 kk=1,6 5 s(kk,j)=s(kk,j)-sp*s(kk,k) sp=spro(j,j,s) 4 call norm(s,j,sp) do 61 ix=1,6 do 61 iy=1,6 61 f6(ix,iy)=f(ind(ix),ind(iy))/dsqrt(m(ix)*m(iy)) c internal force field do 6 ia=1,6 do 6 ib=1,6 fi(ia,ib)=0.0d0 do 6 ix=1,6 do 6 iy=1,6 6 fi(ia,ib)=fi(ia,ib)+f6(ix,iy)*s(ix,ia)*s(iy,ib) fold=fi(1,1) c scaled internal force field do 8 ia=1,6 do 8 ib=1,6 8 fis(ia,ib)=fi(ia,ib)*s6(ia)*s6(ib) fnew=fis(1,1) write(6,602)i,ic(i),io(i),fold,fnew 602 format(i4,':',i4,' - ',i4,', k =',2f12.4) c back transformation: do 7 ix=1,6 do 7 iy=1,6 f6(ix,iy)=0.0d0 do 7 ia=1,6 do 7 ib=1,6 7 f6(ix,iy)=f6(ix,iy)+fis(ia,ib)*s(ix,ia)*s(iy,ib) do 62 ix=1,6 do 62 iy=1,6 62 f(ind(ix),ind(iy))=f6(ix,iy)*dsqrt(m(ix)*m(iy)) 1 continue return end subroutine norm(s,i,sp) implicit none real*8 s(6,6),sp integer*4 i,k do 3 k=1,6 3 s(k,i)=s(k,i)/dsqrt(sp) return end function spro(i,j,s) implicit none integer*4 i,j,k real*8 s(6,6),t,spro t=0.0d0 do 1 k=1,6 1 t=t+s(k,i)*s(k,j) spro=t return end SUBROUTINE READFF(N,FCAR) IMPLICIT none integer*4 N,N1,N3,LN,I,J real*8 FCAR(N,N) OPEN(20,FILE='FILE.FC',STATUS='OLD') N1=1 1 N3=N1+4 IF(N3.GT.N)N3=N DO 130 LN=N1,N 130 READ(20,17)(FCAR(LN,J),J=N1,MIN(LN,N3)) N1=N1+5 IF(N3.LT.N)GOTO 1 17 FORMAT(4X,5D14.6) C DO 31 I=1,N DO 31 J=I+1,N 31 FCAR(I,J)=FCAR(J,I) CLOSE(20) WRITE(*,*)' Cartesian FF read in ... ' RETURN END SUBROUTINE WRITEFF(N,FCAR) IMPLICIT none integer*4 N,N1,N3,LN,J real*8 FCAR(N,N) OPEN(20,FILE='FILES.FC') N1=1 1 N3=N1+4 IF(N3.GT.N)N3=N DO 130 LN=N1,N 130 WRITE(20,17)LN,(FCAR(LN,J),J=N1,MIN(LN,N3)) N1=N1+5 IF(N3.LT.N)GOTO 1 17 FORMAT(I4,5D14.6) CLOSE(20) RETURN END subroutine report(s) character*(*) s write(6,*)s stop end subroutine defa(nat,r,z,ib,ic,io,kb,sf) c define C=O bonds implicit none integer*4 nat,z(*),ib,ic(*),io(*),i,j,kb(*),k,l,nn,no,nc,nh,n real*8 r(*),xi,yi,zi,xj,yj,zj,d,xk,yk,zk,xl,yl,zl,sf(*) character*10 s(5) ib=0 do 1 i=1,nat if(z(i).eq.8)then xi=r(1+3*(i-1)) yi=r(2+3*(i-1)) zi=r(3+3*(i-1)) c Y - X - C = O c l k j i do 2 j=1,nat if(z(j).eq.6)then xj=r(1+3*(j-1)) yj=r(2+3*(j-1)) zj=r(3+3*(j-1)) d=dsqrt((xi-xj)**2+(yi-yj)**2+(zi-zj)**2) if(d.lt.1.34d0)then ib=ib+1 if(ib.gt.3*nat)call report('ib > 3*nat') ic(ib)=j io(ib)=i c atoms bonded to C: c nc .. number of carbons c no .. number of oxygens, except for i c nn .. number of nitrogens c nh .. number of hydrogens on the nitrogen, if any nc=0 nn=0 nh=0 no=0 do 3 k=1,nat if(k.ne.i.and.k.ne.j)then xk=r(1+3*(k-1)) yk=r(2+3*(k-1)) zk=r(3+3*(k-1)) d=dsqrt((xk-xj)**2+(yk-yj)**2+(zk-zj)**2) if(d.lt.1.5d0.and.z(k).eq.7)then nn=nn+1 do 4 l=1,nat if(l.ne.i.and.l.ne.j.and.l.ne.k)then xl=r(1+3*(l-1)) yl=r(2+3*(l-1)) zl=r(3+3*(l-1)) d=dsqrt((xk-xl)**2+(yk-yl)**2+(zk-zl)**2) if(d.lt.1.2d0.and.z(l).eq.1)nh=nh+1 endif 4 continue endif if(d.lt.1.63d0.and.z(k).eq.6)nc=nc+1 if(d.lt.1.63d0.and.z(k).eq.8)no=no+1 endif 3 continue c uknown: kb(ib)=1 c -HN-CO-C: if(nn.eq.1.and.nh.eq.1.and.nc.eq.1)kb(ib)=2 c H2N-CO-C: if(nn.eq.1.and.nh.eq.2.and.nc.eq.1)kb(ib)=3 c =N-CO-C: if(nn.eq.1.and.nh.eq.0.and.nc.eq.1)kb(ib)=4 c O-CO-C: if(no.eq.1.and.nn.eq.0.and.nc.eq.1)kb(ib)=5 endif endif 2 continue endif 1 continue write(6,601)ib 601 format(/,i4,' of C=O bonds:') s(1)='unknown ' s(2)='-NH-CO-C ' s(3)='NH2-CO-C ' s(4)='=N-CO-C ' s(5)='O-CO-C ' do 5 i=1,5 n=0 do 6 j=1,ib 6 if(kb(j).eq.i)n=n+1 write(6,602)n,s(i),sf(i) 602 format(i4,' of ',a10,', sf = ',f10.3) do 7 j=1,ib if(kb(j).eq.i) write(6,603)ic(j),io(j) 603 format(i5,' -',i5,$) 7 if(mod(j,6).eq.0)write(6,*) 5 if(mod(j,6).ne.0)write(6,*) return end