PROGRAM STRETCH
C     stretch molecule according to a normal mode displacement
      IMPLICIT REAL*8 (A-H,O-Z)                                       
      PARAMETER (N0=600)
      REAL*8 SV(N0,3),R(3,N0),Z(N0),FRQ(3*N0),GR(3,N0)
      CHARACTER*80 FILENAME
      CHARACTER*1 OK
      LOGICAL LE
c
      WRITE(*,*)' F.INP - input file'
      WRITE(*,*)
      OPEN(17,FILE='F.INP',STATUS='OLD')                       
      READ(17,*)NQ,NAC
      NOAT=NAC/3
      DO 55 K=1,NOAT                                                  
55    READ(17,*)Z(K),(R(I,K),I=1,3)                                   
      CLOSE(17)
      WRITE(*,*)NOAT,' atoms'
      WRITE(*,6000)1,NQ
6000  FORMAT(' Which modes to stretch (',I1,' - ',I4,'), NMIN,NMAX:')
      READ(*,*)NMIN,NMAX
      WRITE(*,*)'q_max=2 as maximum (y/n)? '
      READ(*,'(a)')OK
      SMAX=2.0D0
      IF(OK.EQ.'n'.OR.OK.EQ.'N')THEN
       WRITE(*,*)'maximum stretch: '
       READ(*,*)SMAX
      ENDIF
      WRITE(*,*)'stretch 0.1 0.05 , -0.5 and -0.1'
      WRITE(*,*)'output in the format of gaussian input'
      WRITE(*,*)'output filename: q1.inp'
      OPEN(3,FILE='q1.inp')
      WRITE(3,1600)
1600  FORMAT('%chk=q1.chk',/,
     1       '#BPW91/6-31G** SP SCF=TIGHT NOSYMM',/,
     2       ' GEOM=(NODISTANCE,NOANGLE)',/,
     2/,
     3'zero geometry',/,
     4/,
     5'0 1')
      DO 52 K=1,NOAT
52    WRITE(3,3000)INT(Z(K)),(R(I,K),I=1,3)
3000  FORMAT(I2,3F17.8)
      WRITE(3,*)
c
      DO 1000 NS=NMIN,NMAX
      DO 1000 ID=1,4
      IF(ID.EQ.1)DELSMALLQ= SMAX
      IF(ID.EQ.2)DELSMALLQ= SMAX*0.5D0
      IF(ID.EQ.3)DELSMALLQ=-SMAX*0.5D0
      IF(ID.EQ.4)DELSMALLQ=-SMAX
c
      OPEN(17,FILE='F.INP',STATUS='OLD')                       
      READ(17,*)NQ,NAC
      DO 35 K=1,NOAT                                                  
35    READ(17,*)Z(K),(R(I,K),I=1,3)                                   
      READ(17,*)                                                      
      DO 45 L=1,NOAT                                                 
      DO 45 I=1,NQ                                                    
      IF(I.EQ.NS)THEN
       READ(17,*)IT,IMD,(SV(L,J),J=1,3)
      ELSE
       READ(17,*)
      ENDIF
45    CONTINUE                                                        
      READ(17,*)                                                      
      READ(17,1717)(FRQ(I),I=1,NQ)                                    
1717  FORMAT(6F11.3)                                                  
      CLOSE(17)
      WS=FRQ(NS)
      DELBIGQ=DELSMALLQ/SQRT(dabs(WS)/219470.0D0)
      DELR=DELBIGQ*0.02342179D0*0.529177249D0
      DO 1 K=1,NOAT
      DO 1 I=1,3
1     R(I,K)=R(I,K)+SV(K,I)*DELR
c
      WRITE(3,1602)NS,WS,DELSMALLQ
1602  FORMAT('--Link1--',/,
     1       '%chk=q1.chk',/,
     2       '#BPW91/6-31G** SP SCF=TIGHT GEOM=(NODISTANCE,NOANGLE)',/,
     3       'NOSYMM GUESS=CHECKPOINT',/,
     4/,
     5'Mode',I4,'(',F8.1,'cm-1), stretched del q = ',F8.3,/,
     6/,
     7'0 1')
      DO 2 K=1,NOAT
2     WRITE(3,3000)INT(Z(K)),(R(I,K),I=1,3)
1000  WRITE(3,*)
c
      CLOSE(3)
      INQUIRE(FILE='FORCES.G94',EXIST=LE)  
      IF(LE)THEN
       WRITE(*,*)' FORCES.G94 exists'
       OPEN(3,FILE='FORCES.G94')
       DO 4 I=1,NOAT
4      READ(3,*)ID,JD,(GR(J,I),J=1,3)
       CLOSE(3)
       GRQ=0.0D0  
       DO 3 I=1,NOAT
       DO 3 J=1,3
3      GRQ=GRQ+GR(J,I)*SV(I,J)
       WRITE(*,6002)GRQ
6002   FORMAT(' Gradient = ',F15.8)
      ENDIF  
      STOP
      END