PROGRAM modesp implicit none character*80 s80 integer*4 nm,n,nat,n1,n2,iargc,i real*8,allocatable::r(:),s(:,:),e(:),m(:),F(:,:),G(:,:), 1p(:,:) integer*4,allocatable::z(:),ml(:) logical list,exclusive write(6,600) 600 format('Projects out selected modes of F.INP from FILE.FC') if(iargc().lt.2)then write(6,6001) 6001 format(' usage: modesp ',/, 1 ' or: modesp <-mmin> (exclusive)',/, 1 ' or: modesp <-m1> ... (list)',/, 1 ' or: modesp <-Nm> <-m1> ... (exclusive list)') stop endif call getarg(1,s80) read(s80,*)n1 exclusive=n1.lt.0 n1=abs(n1) call getarg(2,s80) read(s80,*)n2 list=n2.lt.0 n2=abs(n2) if(list)then allocate(ml(n1)) ml(1)=n2 do 1 i=2,n1 call getarg(1+i,s80) 1 read(s80,*)ml(i) else allocate(ml(n2-n1+1)) endif open(9,file='F.INP') read(9,*)nm,n,nat close(9) allocate(r(n),s(n,n),e(n),m(nat),z(nat),F(n,n),G(n,n),p(n,n)) call reads(n,s,e,nm,'F.INP',r,z) call rm(m,nat) CALL READFF(N,F) CALL mw(N,nm,m,F,G,s,p,0) call PROJECT(G,p,N,n1,n2,e,ml,exclusive,list,nm) CALL mw(N,nm,m,F,G,s,p,1) CALL WRITEFF(N,F) end subroutine mw(N,nm,m,F,G,s,p,ic) implicit none integer*4 N,ic,i,j,nm real*8 m(*),F(N,N),G(N,N),s(N,N),p(N,N),mi if(ic.eq.0)then do 1 i=1,N mi=dsqrt(m((i+2)/3)) do 2 j=1,nm 2 p(i,j)=s(i,j)*mi do 1 j=1,N 1 G(i,j)=F(i,j)/(mi*dsqrt(m((j+2)/3))) else do 3 i=1,N mi=dsqrt(m((i+2)/3)) do 4 j=1,nm 4 s(i,j)=p(i,j)*mi do 3 j=1,N 3 F(i,j)=G(i,j)*mi*dsqrt(m((j+2)/3)) endif return end subroutine rm(m,nat) implicit none real*8 m(*),BM integer*4 nat,NOSC,I OPEN(15,FILE='FTRY.INP') READ(15,*) READ(15,*)NOSC DO 222 I=1,NOSC 222 READ(15,*) READ(15,*) (m(I),I=1,nat) READ(15,*) READ(15,1040)(m(I),I=1,nat) 1040 FORMAT(6G12.6) BM=0.0d0 DO 102 I=1,nat 102 BM=BM+M(I) WRITE(6,778)BM 778 FORMAT(' TOTAL MASS ',F20.6) close(15) return end SUBROUTINE PROJECT(F,s,N,n1,n2,e,ml,exclusive,list,nm) IMPLICIT none real*8 F(N,N),s(N,N),e(*) integer*4 N,n1,n2,ia,ix,i,j,k,n6,nat,ml(*),nm,ip,ic real*8,allocatable::TEM(:,:),A(:,:),A6(:,:) logical list,exclusive integer*4,allocatable::mp(:) nat=N/3 C if(list)then if(exclusive)then WRITE(6,610) 610 format('Conserved modes, E (cm-1):') n6=nm-n1 allocate(mp(n6)) ip=0 do 5 I=1,nm ic=0 do 6 J=1,n1 6 if(ml(J).eq.I)ic=1 if(ic.eq.0)then c this mode is not in list of conserved modes, project it: ip=ip+1 mp(ip)=I endif 5 continue else n6=n1 allocate(mp(n6)) WRITE(6,607) 607 format('Projected modes, E (cm-1):') do 9 i=1,n1 9 mp(i)=ml(i) endif do 1 I=1,n1 1 WRITE(6,609)I,ml(I),e(ml(I)) 609 format(2i5,f10.2) else if(exclusive)then n6=nm-(n2-n1+1) WRITE(6,611)n2-n1+1,e(n1),e(n2),n6 611 format(i4,' modes conserved ',f12.2,' -',f12.2,' cm-1',/, 1 i4,' modes projected') allocate(mp(n6)) ip=0 do 7 I=1,nm c this mode is outside interval of conserved modes, project it: if(I.lt.n1.or.I.gt.n2)then ip=ip+1 mp(ip)=I endif 7 continue else n6=n2-n1+1 WRITE(6,608)N6,e(n1),e(n2) 608 format(i4,' modes projected ',f12.2,' -',f12.2,' cm-1') allocate(mp(n6)) do 8 I=n1,n2 8 mp(I-n1+1)=I endif endif allocate(A6(N,n6),TEM(N,N),A(N,N)) do 21 ia=1,nat do 21 ix=1,3 I=ix+3*(ia-1) do 21 k=1,n6 21 A6(I,k)=s(I,mp(k)) call ORT(A6,n6,nat) DO 22 I=1,N DO 22 J=1,N A(I,J)=0.0d0 DO 22 k=1,N6 22 A(I,J)=A(I,J)-A6(I,k)*A6(J,k) DO 23 I=1,N 23 A(I,I)=A(I,I)+1.0d0 DO 3 I=1,N DO 3 J=1,N TEM(I,J)=0.0d0 DO 3 k=1,N 3 TEM(I,J)=TEM(I,J)+A(I,k)*F(k,J) DO 4 I=1,N DO 4 J=1,N F(I,J)=0.0d0 DO 4 k=1,N 4 F(I,J)=F(I,J)+TEM(I,k)*A(J,k) RETURN END SUBROUTINE ORT(A,N6,NAT) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION A(3*NAT,N6) AMACH=1.0d-13 NAT3=3*NAT C=0.0d0 DO 1 I=1,NAT3 1 C=C+A(I,1)**2 IF(C.GT.AMACH)C=1.0d0/SQRT(C) DO 2 I=1,NAT3 2 A(I,1)=A(I,1)*C DO 3 IC=2,N6 C Column IC, orthonormalize to previous columns: DO 3 ICP=1,IC-1 SP=0.0d0 DO 5 J=1,NAT3 5 SP=SP+A(J,IC)*A(J,ICP) C Subtract projection to ICP: DO 6 J=1,NAT3 6 A(J,IC)=A(J,IC)-SP*A(J,ICP) C Normalize : SP=0.0d0 DO 7 J=1,NAT3 7 SP=SP+A(J,IC)**2 IF(SP.GT.AMACH)SP=1.0d0/SQRT(SP) DO 3 J=1,NAT3 3 A(J,IC)=A(J,IC)*SP 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='FILEP.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 reads(N3,S,E,NQ,fn,r,z) IMPLICIT none integer*4 N3,NQ,nat,i,J,ix,z(*) real*8 S(N3,N3),E(*),r(*) character*(*) fn C 1 a.u.= 2.1947E5 cm^-1 open(4,file=fn) read(4,*)NQ,nat,nat do 1 i=1,NAT 1 read(4,*)z(i),(r(3*(i-1)+ix),ix=1,3) read(4,*) DO 2 I=1,NAT DO 2 J=1,NQ 2 read(4,*)(s(3*(i-1)+ix,NQ-J+1),ix=1,2), 1(s(3*(i-1)+ix,NQ-J+1),ix=1,3) read(4,*) READ(4,4000)(E(NQ-J+1),J=1,NQ) 4000 FORMAT(6F11.6) close(4) write(6,*)NQ,' modes found in '//fn RETURN end