PROGRAM scaleask implicit none integer*4 n,nat,i,nsf,m,j,iold real*8 t real*8,allocatable::F(:,:),B(:,:),sf(:),sf1(:) integer*4,allocatable::ico(:),ico1(:) logical lsep write(6,600) 600 format(/,' Scales selected internal coordinates in', 1' Cartesian force field',/,/, 1' Input: FILE.FC force field',/ 1' B.MAT B-matrix',/, 1' FILE.X geometry',/, 1' ASK.LST list of coordinates and SFs',/,/, 1' Output: FILES.FC scaled force field',/,/) m=1 n=1 allocate(B(m,n)) call rb(0,m,n,B) deallocate(B) allocate(B(m,n)) call rb(1,m,n,B) nat=n/3 write(6,*)'B.MAT read,',nat,' atoms,',m,' coordinates' allocate(ico(m),sf(m),ico1(m),sf1(m)) open(9,file='ASK.LST',status='old') nsf=0 iold=0 sf=1.0d0 lsep=.false. 101 read(9,*,end=99,err=99)t,i if(i.lt.0)then c separate scaling: lsep=.true. i=-i endif do 1 j=iold+1,i nsf=nsf+1 ico(nsf)=j 1 sf(nsf)=t iold=i goto 101 99 close(9) write(6,*)nsf,' scale factors read:' do 2 i=1,nsf 2 write(6,601)i,ico(i),sf(i) 601 format(2i4,f8.4) if(nsf.eq.0)stop allocate(F(N,N)) CALL READFF(N,F) if(lsep)then do 3 i=1,nsf ico1(1)=ico(i) sf1(1)=sf(i) 3 call scale(m,n,1,ico1,sf1,f,B) else call scale(m,n,nsf,ico,sf,f,B) endif CALL WRITEFF(N,F) end subroutine scale(m,n,nsf,ico,sf,f,BB) implicit none integer*4 n,i,j,k,kk,m,nsf,ico(m),a,b,nt real*8 sf(m),f(n,n),sp,spro,BB(m,n),fold real*8,allocatable::T(:,:),fi(:,:),ft(:,:),bt(:,:),bti(:,:), 1fid(:,:) if(nsf.gt.n)call report('nsf>n!') allocate(t(n,n+nsf),bt(nsf,nsf),fid(nsf,nsf),bti(nsf,nsf)) write(6,*)'Old and new force constants:' c define temporary non-redundant B-matrix T, I' = T X c I' ... non-redundant internal coordinates T=0.0d0 c replace part by the desired coordinates, c note that B: Cartesian last, T: Cartesian first do 3 i=1,nsf do 3 j=1,n 3 T(j,i)=BB(ico(i),j) nt=0 do 2 i=1,n 2 T(i,nsf+i)=1.0d0 c orthonormalize until nothing left do 4 j=1,nsf+n do 5 k=1,j-1 sp=spro(j,k,T,n,nsf) do 5 kk=1,n 5 T(kk,j)=T(kk,j)-sp*T(kk,k) sp=spro(j,j,T,n,nsf) if(sp.gt.1.0d-4)then nt=nt+1 call norm(T,j,sp,n,nsf) else do 101 kk=1,n 101 T(kk,j)=0.0d0 endif 4 continue if(nt.ne.n)then write(6,600)nt 600 format(' Rest dimension:',i4) call report('Inconsistent with the number of atoms') endif do 102 j=1,nsf+n if(spro(j,j,T,n,nsf).eq.0)then do 103 k=j,nsf+n-1 do 103 kk=1,n 103 T(kk,k)=T(kk,k+1) endif 102 continue c make non-redundant internal force field F' = T F T allocate(ft(nt,n),fi(nt,nt)) ft=0.0d0 do 6 a=1,nt do 6 j=1,n do 6 i=1,n 6 ft(a,j)=ft(a,j)+t(i,a)*f(i,j) fi=0.0d0 do 7 a=1,nt do 7 b=1,nt do 7 j=1,n 7 fi(a,b)=fi(a,b)+t(j,b)*ft(a,j) c transformation matrix between original (I) and non-redundant (I') c internal coordinates I = BT I' bt=0.0d0 do 8 a=1,nsf do 8 b=1,nsf do 8 i=1,n 8 bt(a,b)=bt(a,b)+BB(ico(a),i)*T(i,b) call INV(nsf,bt,bti) c bti ... laste index redundant c force field in desired redundant coordinates F = BT^-1 F' BT^-1 ft=0.0d0 do 9 a=1,nsf do 9 j=1,nsf do 9 i=1,nsf 9 ft(a,j)=ft(a,j)+bti(i,a)*fi(i,j) fid=0.0d0 do 10 a=1,nsf do 10 b=1,nsf do 10 j=1,nsf 10 fid(a,b)=fid(a,b)+bti(j,b)*ft(a,j) c scaled internal red. force field do 11 a=1,nsf fold=fid(a,a) do 12 b=1,nsf 12 fid(a,b)=fid(a,b)*sf(a)*sf(b) 11 write(6,602)a,fold,fid(a,a) 602 format(i4,2f12.4) c back transformation, to non-red F' = BT F BT: ft=0.0d0 do 13 a=1,nsf do 13 j=1,nsf do 13 i=1,nsf 13 ft(a,j)=ft(a,j)+bt(i,a)*fid(i,j) do 14 a=1,nsf do 14 b=1,nsf fi(a,b)=0.0d0 do 14 j=1,nsf 14 fi(a,b)=fi(a,b)+bt(j,b)*ft(a,j) c back transformation, to Cartesian F = T F' T ft=0.0d0 do 15 a=1,n do 15 j=1,n do 15 i=1,n 15 ft(a,j)=ft(a,j)+t(a,i)*fi(i,j) f=0.0d0 do 16 a=1,n do 16 b=1,n do 16 j=1,n 16 f(a,b)=f(a,b)+t(b,j)*ft(a,j) return end subroutine norm(s,i,sp,n,nsf) implicit none integer*4 i,k,n,nsf real*8 s(n,n+nsf),sp do 3 k=1,n 3 s(k,i)=s(k,i)/dsqrt(sp) return end function spro(i,j,s,n,nsf) implicit none integer*4 i,j,n,k,nsf real*8 s(n,n+nsf),t,spro t=0.0d0 do 1 k=1,n 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 rb(ic,NOB0,NA,B) implicit none integer*4 ic,NA,i,NOB0,nn,ii,ia,ix real*8 b(NOB0,NA) open(9,file='B.MAT') read(9,*)NOB0 read(9,*)na if(ic.eq.0)goto 99 b=0.0d0 read(9,*) do 13 i=1,NOB0 read(9,*) read(9,*)nn do 13 ii=1,nn 13 read(9,*)ia,ix,b(i,ix+3*(ia-1)) 99 close(9) return end SUBROUTINE INV(n,a,ai) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) dimension ai(n,n),a(n,n) real*8,allocatable::e(:,:) C allocate(e(n,2*n)) TOL=1.0d-9 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 write(6,*)' line ',ii write(6,188)tol 188 format(' tol = ',e12.4) write(6,*)'Cannot invert matrix' stop 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