PROGRAM a162cgGS5 C CCC for poly ALA15 potential N=486 CCCCCCCCCCCCCCCCCCCCCCCCCC C W.Quapp 30.11.2004 /07.07.2006 for RGF=Newton Trajectory C GS nach: Baron Peters et al. JCP 120 (2004) 7877-7886 c script version - modular structure CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC c predictor step with z-Mat coordinates CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC Integer N,L,N3 PARAMETER(N=486,L=12,N3=162,N6=480) REAL*8 Y(L,N),Ye(N),Yint(L,N6),ENall(L) Integer J,I,LA,ITall,Istatus,itGS,iflag,Natoms(N3) LOGICAL FINISH Character*2 CharAtoms(N3) Character*50 zmatLinesLA1(N3+4),zmatLinesYe(N3+4) Character*7 znames(N6) Data CharAtoms/ 1 'C','C','N','C','C','O','C','N','C','C','O','C','N','C','C','O', 1 'C','N','C','C','O','C','N','C','C','O','C','N','C','C','O','C', 1 'N','C','C','O','C','N','C','C','O','C','N','C','C','O','C','N', 1 'C','C','O','C','N','C','C','O','C','N','C','C','O','C','N','C', 1 'C','O','C','N','C','C','O','C','N','C','C','O','C','N','C','O', 1 'H','H','H','H','H','H','H','H','H','H','H','H','H','H','H','H', 1 'H','H','H','H','H','H','H','H','H','H','H','H','H','H','H','H', 1 'H','H','H','H','H','H','H','H','H','H','H','H','H','H','H','H', 1 'H','H','H','H','H','H','H','H','H','H','H','H','H','H','H','H', 1 'H','H','H','H','H','H','H','H','H','H','H','H','H','H','H','H', 1 'H','H'/ c Data Natoms/6,6,7,6,6,8,6,7,6,6,8,6,7,6,6,8, * 6,7,6,6,8,6,7,6,6,8,6,7,6,6,8,6, * 7,6,6,8,6,7,6,6,8,6,7,6,6,8,6,7, * 6,6,8,6,7,6,6,8,6,7,6,6,8,6,7,6, * 6,8,6,7,6,6,8,6,7,6,6,8,6,7,6,8, * 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, * 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, * 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, * 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, * 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, * 1,1/ cc c 7: all points of chain for Mma input OPEN(7,FILE='a162gs.weg') c 77: all energy points for Mma input OPEN(77,FILE='a162keten.weg') c 9: all points of chain OPEN(9,FILE='a162chain.dat') OPEN(11,FILE='a162param.dat') OPEN(12,FILE='a162proj.dat') c 13: actual point X of GS iteration OPEN(13,FILE='a162point.dat') OPEN(14,FILE='a162ener.dat') c OPEN(15,FILE='a162grad.dat') OPEN(21,FILE='a162iflag.dat') OPEN(22,FILE='a162cgplFi.dat') OPEN(23,FILE='a162cggsW.dat') cc OPEN(32,FILE='a162atoms.dat',status='old') OPEN(33,FILE='a162poig.dat') OPEN(44,FILE='a162proto.txt',status='unknown',access='append') cc z-mat coordinates OPEN(83,FILE='starta162.zmat') OPEN(84,FILE='endea162.zmat') OPEN(85,FILE='a162current.zmat') OPEN(89,FILE='a162chainint.zmat') cccc c GS has actual point at node LA on chain to final minimum C Constants c N=66 c N3=N/3 c N6=n-6 IFLAG=0 ww=0.0d0 rewind 21 write(21,*) IFLAG rewind 23 write(23,*) (ww,I=1,N+1) rewind 22 FINISH=.FALSE. write(22,*) FINISH, IFLAG, IFLAG istatus=0 c itGS iteration counter in the GS loop rewind 11 READ(11,*) J, LA, EPS, ITALL, itGS J=L WRITE(44,*) L, LA, EPS, ITALL, itGS WRITE(44,*) ' predictor step to LA+1 = ',LA+1 c c Notbremse: emergency brake only c IF(itGS.gt. 3*N+3) istatus=10 c rewind 9 Do 11, J=1,L READ(9,*) J1 READ(9,*) Do 11, JJ=1,N3 11 READ(9,555) CharAtoms(JJ),(Y(J,3*(JJ-1)+I),I=1,3) 555 Format(1X,A3, 3F16.8) c rewind 14 READ(14,*) ENall(LA) WRITE(44,*) ' energy before predictor ', ENall(LA) IF( LA.eq.1) WRITE(77,*) ENall(LA) If( LA.gt.1) then rewind 77 READ(77,*) (ENall(I),I=1,LA-1) rewind 77 WRITE(77,*)(ENall(I),I=1,LA) ENDIF rewind 13 READ(13,*)(Y(LA,I),I=1,N) c c GS searches actual point at node LA+1 on chain to final minimum ccccccccccccccccccccccccccccccccccccccccccccc c ! c LA=LA+1 c ! c ccccccccccccccccccccccccccccccccccccccccccccc rewind 11 WRITE(11,*) L, LA, EPS ,ITALL , itGS WRITE(44,*) L, LA, EPS, ITALL, itGS WRITE(44,*) ' Predictor: go with LA one step further ' ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c new guess point in Cartesian coordinates C extra besides the script organization: Do 712, J=1,L-LA+1 Do 713, I=1,N 713 Y(LA-1+J,I)=( Y(LA-1,I)*(L-LA-J+1)+ Y(L,I)*J )/(L-LA+1) 712 continue rewind 9 Do 714, J=1,L WRITE(9,*) ' ',N3 WRITE(9,*) ' ' Do 714, JJ=1,N3 714 WRITE(9,555) CharAtoms(JJ),(Y(J,3*(JJ-1)+I),I=1,3) rewind 13 write(13,50)(Y(LA,I),I=1,N) rewind 33 Do 17, J=1,N3 17 WRITE(33,*) Natoms(J), (Y(LA,3*(J-1)+I),I=1,3) WRITE(33,*) ' ' ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc If(LA.eq.L) then rewind 13 write(13,50)(Y(L,I),I=1,N) ITall=ITall+itGS WRITE(6,25) LA, ITall itGS=0 rewind 11 WRITE(11,*) L, LA, EPS ,ITALL , itGS istatus=10 WRITE(6, *) ' --- STOP: Final Minimum --- ' c Mma - output of Reaction Path rewind 7 WRITE(7,*) '{' Do 15, J=1,L-1 WRITE(7,22)(Y(J,I),I=1,N3) 15 continue WRITE(7,21)(Y(L,I),I=1,N3) rewind 77 READ(77,*) (ENall(I),I=1,L) rewind 77 WRITE(77,*) ' { ' WRITE(77,23)(ENall(I),I=1,L-1) WRITE(77,*) ENall(L), ' } ' WRITE(6, *) ' --- Sum of Gradient Calculations -- >> ',ITall WRITE(44,*) ' --- Sum of Gradient Calculations -- >> ',ITall goto 999 ENDIF c rewind 85 Do 3 I=1,N3+1 3 READ(85,833) zmatLinesLA1(I) 833 Format(50A) Do 333 I=1,N6 333 READ(85,*) znames(I),Yint(LA-1,I) Do 4 I=1,N3+1 READ(84,833) zmatLinesYe(I) cccc control order of dist, angles, dihedrals ccccccccccccccccc If(zmatLinesLA1(I). ne. zmatLinesYe(I)) then write(6,*) 'ERROR in z-mat structure is to correct by hand !!' write(6,*) ' LINE = ', I write(44,*)'ERROR in z-mat structure is to correct by hand !!' write(44,*)' LINE = ', I endif 4 continue Do 41 I=1,N6 41 READ(84,*) znames(I),Ye(I) c harmonize the sign of the dihedrals: Do 311 I=6,N6,3 IF(Yint(LA-1,I)*Ye(I) .lt. 0.0d0 1 .and. Dabs(Yint(LA-1,I)).gt. 130.0d0 2 .and. Dabs( Ye(I)).gt. 130.0d0 ) then write(44,*)'dih - ERROR in Yint(',LA-1,I,') corrected ' if (Yint(LA-1,I).gt.0.0d0 ) then Yint(LA-1,I)=Yint(LA-1,I)-360.0d0 else Yint(LA-1,I)=Yint(LA-1,I)+360.0d0 endif ENDIF 311 continue c c new guess points Do 12 J=1,L-La Do 111 I=1,N6 111 Yint(LA-1+J,I)= ((L-LA-J+1)*Yint(LA-1,I)+J*Ye(I))/(L-LA+1) 12 continue rewind 89 Do 2, J=LA,L-1 Do 32 I=1,N3+1 32 Write(89,833) zmatLinesLA1(I) Do 321 I=1,N6 321 WRITE(89,*) znames(I), Yint(J,I) WRITE(89,*) 2 continue Do 33 I=1,N3+1 33 Write(89,833) zmatLinesYe(I) Do 322 I=1,N6 322 WRITE(89,*) znames(I), Ye(I) WRITE(89,*) ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc rewind 85 Do 3332 I=1,N3+1 3332 Write(85,833) zmatLinesLA1(I) Do 3321 I=1,N6 3321 WRITE(85,*) znames(I), Yint(LA,I) WRITE(85,*) ITall=ITall+itGS itGS=0 WRITE(6,25)LA,ITall WRITE(44,25)LA,ITall rewind 11 WRITE(11,*) L, LA, EPS ,ITALL , itGS WRITE(44,*) L, LA, EPS ,ITALL , itGS C end iteration cycle 999 CONTINUE 21 FORMAT(3H { , 21(F15.9, 3H , ), F15.9, 3H }} ) 22 FORMAT(3H { , 21(F15.9, 3H , ), F15.9, 3H }, ) 23 FORMAT(F15.9, 3H , ) 25 FORMAT(' Node ',I4,' ITall ',I5) 50 FORMAT(1X,3F20.10) close(44) call exit(istatus) STOP END