PROGRAM A22ntGS3 C CCC for Ala potential NDIM=60 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C W.Quapp June 2007 for Newton Trajectory by GS following C GS nach: Baron Peters et al. JCP 120 (2004) 7877-7886 C W.Q.: JCP 122, iss17, (2005) 174106 (11 pages) C c script version for GamessUS (or other PES function) c GS cycle for RP: NewPoint calculation for ProjectedGradient=0 c corrector step is done by second order methods of NT theory c program parts from Michael Hirsch (Dissertation 2003, p.81) C C the old NT point for predictor step is Xold C the actual point, to start the correction, is Y(LA,..) C put tangent = X-Yold = X - Y(LA-1) Ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc Integer NDIM,L,N3 PARAMETER(NDIM=58,L=30,N3=22) c steplength-suche mit 2308 Integer LA,N,Nm1 Integer istatus,ITALL,Natoms(N3) REAL*8 energy,gnorm,tnorm,EPS,vec_ddot,D1 c REAL*8 torad,btoa REAL*8 VV(NDIM),told(NDIM),X(NDIM),Proj(NDIM) REAL*8 Y(L+1,NDIM),Yold(NDIM),Xold(NDIM) REAL*8 Rs(NDIM),ProMat(NDIM,NDIM),hessian(NDIM,NDIM) c ,gold(NDIM) ,hold(NDIM,NDIM) Integer JJ,I,J,Np,Np1,IterMax,IFLAG,NFUN Character*2 CharAtoms(N3) Character*5 Zcoor(NDIM+2) c CHASS: Data Zcoor/'ch2 ','cc3 ','cch3 ','oc4 ','occ4 ','dih4 ', 1 'hc5 ','hcc5 ','dih5 ','hc6 ','hcc6 ','dih6 ','nc7 ', 2 'ncc7 ','dih7 ','cn8 ','cnc8 ','dih8 ','cc9 ','ccn9 ', 3 'dih9 ','oc10 ','occ10','dih10','hc11 ','hcn11','dih11', 7 'hn12 ','hnc12','dih12','cc13 ','ccn13','dih13','hc14 ', 8 'hcc14','dih14','hc15 ','hcc15','dih15','hc16 ','hcc16', 7 'dih16','nc17 ','ncc17','dih17','cn18 ','cnc18','dih18', 8 'hc19 ','hcn19','dih19','hn20 ','hnc20','dih20','hc21 ', 7 'hcn21','dih21','hc22 ','hcn22','dih22'/ Data CharAtoms/'h','c','c','o','h','h','n','c','c','o','h','h', 1 'c','h','h','h','n','c','h','h','h','h'/ Data Natoms/1,6,6,8,1,1,7,6,6,8,1,1,6,1,1,1,7,6,1,1,1,1/ c OPEN(9,FILE='ala2chain.dat') OPEN(10,FILE='ala2Rs.dat') OPEN(11,FILE='ala2param.dat') OPEN(12,FILE='ala2proj.dat') OPEN(13,FILE='ala2point.dat') OPEN(14,FILE='ala2ener.dat') OPEN(15,FILE='ala2grad.dat') OPEN(16,FILE='ala2engs.weg') OPEN(18,FILE='ala2gmatr.dat') OPEN(19,FILE='ala2prest.dat') OPEN(20,FILE='ala2hesse.dat') OPEN(21,FILE='ala2iflag.dat') OPEN(25,FILE='ala2told.dat') OPEN(27,FILE='ala2xold.dat') OPEN(28,FILE='ala2yold.dat') OPEN(33,FILE='ala2poig.dat') c OPEN(34,FILE='ala2gold.dat') c OPEN(35,FILE='ala2hold.dat') OPEN(44,FILE='ala2protocol.txt',status='unknown',access='append') OPEN(54,FILE='ala2dih917p.txt',status='unknown',access='append') OPEN(55,FILE='ala2dih917a.txt',status='unknown',access='append') c Read problem input information N = NDIM Nm1= NDIM -1 C N3 = (N+6)/3 C N3: number of atoms ITERMAX= 10 istatus= 1 c (bohr,radian)->(angstroem,degree) for internal coordinates c torad=dacos(-1.0d0)/180.0d0 c btoa =0.52917706d0 c read(21,*) IFLAG write( 6,*) ' IFLAG is ', IFLAG write(44,*) ' IFLAG is ', IFLAG c c IFLAG=0 indicates an initial entry into corrector c rewind 11 read(11,*) JJ, LA, eps ,ITALL, NFUN cc JJ=L is constant !! rewind 11 write( *,*) L, LA, eps ,ITALL, NFUN write(11,*) L, LA, eps ,ITALL, NFUN write(44,*)' L=',L,' LA= ',LA,' ITall = ',ITALL,'+',NFUN IF(LA.eq.L+1) then istatus=10 goto 999 Endif c Print parameters if ( LA.lt.3 .AND. IFLAG.eq.0)then write(*,820) write(*,840) N write(44,820) write(44,840) N endif c current point X rewind 13 READ(13,*) (X(I),I=1,N) write(44,*) ' X is ' write(44,50) (X(I),I=1,5) write(44,50) (X(I),I=6,N) c string Y rewind 9 Do 4, J=1,L+1 Read(9,50) (Y(J,I),I=1,N) 4 Read(9,50) c current string Y goes up to LA if(IFLAG.eq.0) then do 45 I=1,N 45 Yold(I) = Y(LA-1,I) rewind 28 write(28,50) (Yold(I),I=1,N) c write(*,*) ' Yold is ' c write(*,50) (Yold(I),I=1,N) Do 211 I=1,N 211 Xold(I)= Y(LA,I) write(44,*) 'IFLAG is ', IFLAG rewind 27 write(27,50) (Xold(I),I=1,N) else read(27,50) (Xold(J),J=1,N) read(28,50) (Yold(J),J=1,N) c read(34,50) (gold(J),J=1,N) c write(44,*) ' Xold is ' c write(44,50) (Xold(I),I=1,N) endif c IF((LA.lt.3) .and. (IFLAG.eq.0)) then do 451 I=1,N 451 told(I) = X(I)- Yold(I) c (bohr,radian)<<---(angstroem,degree) for internal coordinates call angtobohr(told,N) gnorm=DSQRT(vec_ddot(told,told,N)) call vec_dmult_constant(told,N,1.0d0/gnorm,told) write(44,*) ' first told is ' write(44,50) (told(I),I=1,5) write(44,50) (told(I),I=6,N) else read(25,50) (told(J),J=1,N) write(44,*) ' told is read ' write(44,50) (told(I),I=1,5) write(44,50) (told(I),I=6,N) endif cc If(NFUN .gt. Itermax) then write(6,*) ' Cancel Loop: ERROR , no Convergence at LA',LA write(6,*) ' LA= ' ,LA, ' NFUN = ', NFUN write(44,*)' Cancel Loop: ERROR , no Convergence at LA ',LA istatus=1 IFLAG=2 ITALL = ITALL+NFUN NFUN=0 write(11,*) L, LA, eps ,ITALL, NFUN write(44,*)' LA= ',LA,' ITall= ',ITALL,' NFUN= ',NFUN write(44,*) write(44,*) ' STOP corrector ' goto 111 Endif c read the SearchDirection Rs and (old) global ProjectionMatrix here rewind 10 READ(10,50) (Rs(I),I=1,N) if(IFLAG.eq.0) WRITE(44,*) 'Rs ' if(IFLAG.eq.0) WRITE(44,50) (Rs(I),I=1,5) if(IFLAG.eq.0) WRITE(44,50) (Rs(I),I=6,N) gnorm=vec_ddot(Rs,told,N) WRITE(44,*)' Cartesian Rs * told =', gnorm c rewind 12 Do 3 J=1,N READ(12,50) (ProMat(J,I),I=1,N) cc if(IFLAG.eq.0) WRITE(44,*) ' J= ',J cc if(IFLAG.eq.0) WRITE(44,50) (ProMat(J,I),I=1,N) 3 READ(12,50) c read function, chain, and gradient values here rewind 14 READ(14,*) energy WRITE(44,*)'energy=',energy WRITE(*,*) 'energy=',energy c file with Grad rewind 15 READ(15,*) (VV(I),I=1,6) c note: override first dih4 of first CH3 group troughout !! READ(15,*) (VV(I),I=6,N) write(44,*) ' gradient is ' write(44,50) (VV(I),I=1,5) write(44,50) (VV(I),I=6,N) c cccc IF( NFUN .eq.0 ) then rewind 20 Do 120, J=1,6 READ(20,52) ( hessian(J,I),I=1,6 ) 120 READ(20,52) ( hessian(J,I),I=6,N ) Do 121, J=6,N READ(20,52) ( hessian(J,I),I=1,6 ) 121 READ(20,52) ( hessian(J,I),I=6,N ) c test write(44,*) ' Diag of Hessian is ' write(44,50) (hessian(I,I),I=1,5) write(44,50) (hessian(I,I),I=6,N) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Hess Update c call vec_dcopy(hessian,hold,N*N) c write(44,*) ' === start analytical Hessian === ' c cccc ENDIF c IF( NFUN.gt.0 ) then c rewind 35 c Do 122, J=1,N c122 READ(35,50) (hold(J,I),I=1,N) c write(44,*) ' === use updated Hessian from former LA === ' c ENDIF c write(44,*) ' NFUN = ', NFUN c write(44,*) ' hessian matrix befor Corr' c Do 7770, J=1,N c7770 WRITE(44,50) (hold(J,I),I=1,N) c c Cartesian norm of gradient gnorm= DSQRT(vec_DDOT(VV,VV,N)) WRITE(*,*) ' PES Cartesian grad norm ',gnorm WRITE(44,*)' PES Cartesian grad norm ',gnorm ccccccccccccccccccccccccccccccccccccccccccALT: C Newton trajectory with Projector to Rs C Old: do the "Cartesian" Projektion of Gradient call vec_dcopy(VV,Proj,N) call PROJECTION(N,Proj,ProMat) D1= DSQRT(vec_DDOT(Proj,Proj,N)) c WRITE(6,*) ' projected grad norm ', D1 WRITE(44,*)' projected Cart. grad norm', D1 If(D1 .lt. 1.0d0/30.0d0*eps) THEN WRITE(44,*)'at IteGS=',NFUN,' ok Cart proj norm less then', 1 ' 1/30 * tolerance EPS ', eps/30.0d0 WRITE(44,*)' successful convergence (no ERROR) ' WRITE(6, *)' IteGS=',NFUN,' == ok - Cart proj norm',D1 ITALL = ITALL+NFUN NFUN=0 IFLAG=2 istatus=1 GOTO 111 ENDIF c c Call the new Newton corrector code cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc CALL newtcorr(NDIM,IFLAG,LA,Eps,X,VV,hessian,Rs,Xold,told) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c IFLAG= c 1 : Corrector step inside the corr - loop c 2 : return to predictor cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc NFUN=NFUN+1 rewind 11 write(11,*) L, LA, eps ,ITALL, NFUN c actual point approximation: rewind 13 WRITE(13,50) (X(I),I=1,N) ccc for gamessUS input file rewind 33 Do 7, J=1,5 7 WRITE(33,*) Zcoor(J), X(J) Do J=6,N WRITE(33,*) Zcoor(J+1), X(J) enddo WRITE(33,*) ' ' rewind 21 write(21,*) IFLAG c write(44,*)' IFLAG after corr ', IFLAG write(44,*) rewind 25 write(25,50) (told(J),J=1,N) c c call vec_dcopy(VV,gold,N) c rewind 35 c Do 77, J=1,N c77 WRITE(35,50) (hold(J,I),I=1,N) c write(44,*) ' hessian-matrix nach Corr' c Do 770, J=1,N c770 WRITE(44,50) (hold(J,I),I=1,N) c Do 771, J=1,N c771 WRITE(20,50) (hold(J,I),I=1,N) c WRITE(20,50) c rewind 34 c write(34,50) (gold(J),J=1,N) ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc write (*,*) ' energy ', energy IF(IFLAG.EQ.2) THEN rewind 27 write(27,50) (Xold(J),J=1,N) write(44,890) energy write(44,*) Do 17, J=1,5 17 WRITE(44,*) Zcoor(J),' ',X(J) Do J=6,N WRITE(44,*) Zcoor(J+1),' ',X(J) enddo c output for Mma Figure: Do J=1,N If(Zcoor(J).eq.'dih9') WRITE(54,*) '{ ',X(J-1),' ,' If(Zcoor(J).eq.'dih17') WRITE(54,*) X(J-1),' }, ' If(Zcoor(J).eq.'dih9') WRITE(55,*) '{ ',X(J-1),' ,' If(Zcoor(J).eq.'dih17') WRITE(55,*) X(J-1),' }, ' Enddo IF(LA.le.L) write(44,*)' new step of NewtonTrajectory: ' ENDIF ccccccccccccccc 111 continue ccccccccccccccc rewind 11 write(11,*) L, LA, eps ,ITALL, NFUN rewind 21 write(21,*) IFLAG c Istatus controls the corrector loop IF(IFLAG.EQ.2) THEN istatus=1 ELSE istatus=0 ENDIF c new chain point becomes: Do 117 I=1,N 117 Y(LA,I)= X(I) rewind 9 Do 118, J=1,L+1 write(9,50) (Y(J,I),I=1,N) 118 write(9,50) WRITE(16,*) '{ ' IF(LA.le.L) WRITE(16,*) ' ', energy,' , ' IF(LA.eq.L+1)WRITE(16,*) ' ', energy,' } ' 999 continue close(44) call exit(istatus) STOP 50 FORMAT(1X,3F20.10) 51 FORMAT(1X, F25.16) 52 FORMAT(1X,3F15.10) 820 format ( ' *** NT Routine for ProjGrad=zero *** ') 840 format ( ' Ndim =',i6 ) 890 format (/,' Energy at pes(x*) = ', 1pd16.8) END cccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE PROJECTION(N,VEC,ProMat) c projection of vector VEC of dim N c by Projection Matrix ProMat of dim(N,N) c ! old version in Cartesian metric ! REAL*8 VEC(N),ProMat(N,N),Proj(N) INTEGER N,J,JJ Do 1 J=1,N Proj(J)=0.0d0 Do 2 JJ=1,N 2 Proj(J)=Proj(J) + Vec(JJ)*ProMat(J,JJ) 1 Continue Do 3 J=1,N 3 Vec(J) = Proj(J) RETURN END Ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc C subroutine newtcorr(nn,IFLAG,LA,EPS,point,grad,hessian,rsel, * xold,told) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC W.Quapp, February 2008 CC for GS-NT corrector by NT second order method CC with parts of a program written by Michael Hirsch CC use metric of internal z matrix coordinates CC CC Pred: replace tangent of NT methods by predictor step of GS CC ( next program part ***5.f ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC integer i,j,nn,n3 real*8 EPS,EPSneu,hessian,gradcont, c torad,btoa,dsqrt,gold, 1 rsel,dk,gmat,unit,rgeq,rgrad,point,vec_mdot, 2 grad,tang,tcov,told,gnorm,dnst,tgt,tg,dmnorm real*8 xold,projNm1,tt,vec_ddot,dstep,rgnorm,ginv,pla, 3 prestep dimension rsel(nn),dk(nn,nn+1),hessian(nn,nn),grad(nn), 1 rgeq(nn-1,nn),dstep(nn),projNm1(nn-1,nn), 2 rgrad(nn-1),point(nn),xold(nn),gradcont(nn), c gold(nn),hold(nn,nn),dx(nn),dg(nn), 3 tcov(nn),tang(nn),told(nn),gmat(nn,nn),unit(nn,nn), 4 ginv(nn,nn),prestep(nn) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c write(44,*) '( point ) in angstroem/degree in Corrector ' c write(44,50) point(1) c write(44,50) point(2), point(3) c write(44,50) (point(I),I=4,NN) n3=(nn+ 6 +2 )/3 c test holds only for full dimensional case: c tt=(point(2)**2)*Dsin(point(3)) c Do 327 K=2,NAT/3 c tt=tt*(point(3*(K-1)+1)**2)*Dsin(point(3*(K-1)+2)) c 327 continue c write(44,*) 'theor.Sqrt of |g_ij| in angstroem ',tt c tt=tt*(1.0d0/0.52917706d0)**21 c write(44,*) 'theor.Sqrt of |g_ij| in bohr ',tt cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c (bohr,radian)<<--(angstroem,degree) for internal coordinates c note: angtobohr depends from the z-matrix used !! c it is to adapt at every change !! c torad=dacos(-1.0d0)/180.0d0 c btoa =0.52917706d0 call angtobohr(point,NN) write(44,*) '( point ) in bohr,radian ' write(44,50) (point(I),I=1,5) write(44,50) (point(I),I=6,NN) write(44,*) call mat_diag(unit,nn,1.0d0) c c note: gmat and ginv interchanged agains use in Vibrational Theory c here: mathematical gmat is metric g_ij, c g^ij is its contravariate inverse c now gmat is obtained by svd of ginv, ginv is from the B-matrix c rewind 18 Do J=1,6 READ(18,50) (ginv(J,I),I=1,6) READ(18,50) (ginv(J,I),I=7,nn) enddo Do J=7,nn READ(18,50) (ginv(J,I),I=1,6) READ(18,50) (ginv(J,I),I=7,nn) enddo c READ(18,*) Do 767, J=1,6 READ(18,50) (gmat(J,I),I=1,6) 767 READ(18,50) (gmat(J,I),I=7,nn) Do J=7,nn READ(18,50) (gmat(J,I),I=1,6) READ(18,50) (gmat(J,I),I=7,nn) enddo cccccccc write(44,*) ' ginv I,I is ' WRITE(44,50) (ginv(I,I),I=1,5) WRITE(44,50) (ginv(I,I),I=6,nn) write(44,*) c write(44,*) ' gmat(21,I) is dih 9=phi ' c WRITE(44,50) (gmat(21,I),I=1,nn) c write(44,*) ' ginv(45,I) is dih 17=psi' c WRITE(44,50) (ginv(45,I),I=1,nn) c write(44,*) write(44,*) ' gmat I,I is ' WRITE(44,50) (gmat(I,I),I=1,5) WRITE(44,50) (gmat(I,I),I=6,nn) c unit test call matmult(ginv,nn,gmat,nn,unit,nn) write(44,*) ' ginv*gmat I,I is ' WRITE(44,50) (unit(I,I),I=1,nn) write(44,*) ' ginv*gmat 17,I is ' WRITE(44,50) (unit(17,I),I=1,nn) gnorm=dmnorm(grad,ginv,nn) WRITE(44,*) 'PES full grad-metric norm',gnorm ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cc test: it comes out successful c write(44,*)' only DIAG METRIC is used in ginv' c Do 7661, J=1,nn c Do 7661, I=1,nn c7661 If (I.ne.J) ginv(J,I)=0.0d0 c call mat_diag(gmat,nn,1.0d0) c write(44,*)'DIAG METRIC 1/diag(ginv) is used in gmat' c Do 7671, J=1,nn c7671 gmat(J,J)=1.0d0/ginv(J,J) ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc C contravariant gradient call vec_dinit(gradcont,nn,0.0d0) call matmult(ginv,nn,grad,1,gradcont,nn) c write(44,*) ' contra gradient is ' c WRITE(44,50) (gradcont(I),I=1,nn) ccccccccccccccccccccccccccccccccccccccccccccccccccccc c Version of old rgf: direction of rsel with first gradient... c if (istep.eq.0 .and. .... ) then c call matmult(ginv,nn,gradient,1,rsel,nn) c call ortcompl(rsel,proj,unit,nn) c endif c Note: rsel not used in metric form, up to now !! tt=dmnorm(rsel,gmat,nn) WRITE(44,*)'metric norm of rsel befor normalization ',tt call vec_dmult_constant(rsel,nn,1.0d0/tt,rsel) c C rsel search direction to (N-1)xN - ProjMat cc test: gmat to make projNm1 lines contravariate cc call ortcompl(rsel,projNm1,gmat,nn) ccc but use again old version by MH: call ortcompl(rsel,projNm1,unit,nn) c write(44,*)' projNm1 is ' c do 121 J=1,nn-1 c write(44,50) (projNm1(J,I),I=1,nn) c121 write(44,50) c then here unit is correct for Pr * H cc call tri_mat(projNm1,nn-1,nn,unit,nn,hessian,nn,rgeq) ccc but use again old version by MH: call tri_mat(projNm1,nn-1,nn,ginv,nn,hessian,nn,rgeq) write(44,*)' rgeq(I,I) is ' write(44,50) (rgeq(I,I),I=1,5) write(44,50) (rgeq(I,I),I=6,nn-1) call vec_dinit(rgrad,nn-1,0.0d0) cc projNm1 is contravariate, thus it is to apply to grad c gradcont is useless now cc call matmult(projNm1,nn-1,grad,1,rgrad,nn) ccc but use again old version by MH: call matmult(projNm1,nn-1,gradcont,1,rgrad,nn) IF(IFLAG.eq.0) then write(44,*)' rgrad is ' write(44,50) (rgrad(I),I=1,nn-1) ENDIF rgnorm=dsqrt(vec_ddot(rgrad,rgrad,nn-1)) WRITE(6,*) ' projected rgnorm ',rgnorm WRITE(44,*)' projected rgnorm ',rgnorm c Tangent : gmat is used for normalization of the tangent c which is contravariate call gettang(rgeq,gmat,tang,nn) write(44,*) ' tang is ' write(44,50) (tang(I),I=1,5) write(44,50) (tang(I),I=6,nn) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccctest c new (calculation without bifurcation test) tgt=vec_mdot(tang,gmat,told,nn) write(44,*) ' tang*gmat*told is ', tgt if (tgt.lt.0.0d0) call vec_dmult_constant(tang,nn,-1.0d0,tang) call vec_dcopy(tang,told,nn) call vec_dinit(tcov,nn,0.0d0) call matmult(gmat,nn,tang,1,tcov,nn) write(44,*) ' tcov is ' write(44,50) (tcov(I),I=1,5) write(44,50) (tcov(I),I=6,nn) pla=0.0d0 cccccccccccccccccccccccccccccccccccccccccccccc c call buildk(dk,rgeq,tang,rgrad,nn,pla) ccc but use again old version by MH: call buildk(dk,rgeq,tcov,rgrad,nn,pla) cccccccccccccccccccccccccccccccccccccccccccccc c t1 for bifurcation detection in old program: c detk=det(dk,nn) c write(44,*) ' detK is ', detk c if (detk.lt.0.0d0) then c call vec_dmult_constant(tang,nn,-1.0d0,t1) c else c call vec_dcopy(tang,t1,nn) c endif ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c tgt=vec_mdot(tang,gmat,told,nn) c write(44,*) ' tang*gmat*told is ', tgt c tgt=vec_ddot(tang,told,nn) c write(44,*) ' tang 1 * told is ', tgt cc if (tgt.lt.0.0d0) then cc call vec_dmult_constant(tang,nn,-1.0d0,tang) cc call vec_dmult_constant(tcov,nn,-1.0d0,tcov) cc call buildk(dk,rgeq,tcov,rgrad,nn,pla) cc endif cc call vec_dcopy(tang,told,nn) ccccccccccccccccccccccccccccccccccccccccccccc c tighter near stationary points the eps: c EPSneu=EPS c If(Eps.gt.0.1d0*gnorm .and. LA.gt.5) then c EPSneu=0.1d0* gnorm c write(44,*) 'new eps is 0.1*gnorm',EPSneu c else c write(44,*) ' eps is ', EPS c Endif cccccccccccccccccccccccccccccccccccccccccccccc if (rgnorm.gt.EPS) then c Corrector call linsolve(dk,dstep,nn) tt=dmnorm(dstep,gmat,nn) write(*,*) 'test dstep norm tdn ',tt if(tt .gt. 1.0d0) tt=1.0d0 read(19,50) (prestep(I),I=1,nn) call angtobohr(prestep,nn) pla=vec_ddot(prestep,tcov,nn) write(*,*) 'prestep*tcov ',pla tgt=dmnorm(prestep,gmat,nn) write(*,*) ' |predstep| ',tgt pla=pla/tgt tgt=dmnorm(tcov,ginv,nn) write(*,*) ' norm tcov ',tgt pla=pla/tgt write(*,*) 'pla vor delta',pla pla=Dsqrt(Dabs(1.0d0-pla**2))*tt*0.90d0 write(*,*) ' sin(alpha)*tdn*0.9 = ',pla c If(pla .gt. 0.25d0) pla=0.25d0 write(44,*) ' delta von pred-step ', pla call buildk(dk,rgeq,tcov,rgrad,nn,pla) call linsolve(dk,dstep,nn) ccccccccccccccccccccccccccccccccccccccccccccccccccccc IFLAG=1 goto 111 else c goto Predictor IFLAG=2 write(44,*) write(44,*)'Convergence (no ERROR) - go to next Pred, LA=',LA c (bohr,radian)-->>(angstroem,degree) call bohrtoang(point,nn) call vec_dcopy(point,xold,nn) endif return ccccccccccccccccccccccccccccccccccccccccccccccccccc c normalization of correctors dstep 111 continue tt=dmnorm(dstep,gmat,nn) ttr=tt write(44,*) ' step-norm of corrector ', tt stepmax=1.000000000000000000002308d0 if(tt .lt. stepmax) then tt=1.0d0 else write(44,*)' ! ERROR corr-step shorten to ',stepmax tt=stepmax/tt endif c execute step call vec_dmult_add(point,dstep,nn,tt,point) write(44,*) ' dstep is ' write(44,50) (tt*dstep(I),I=1,5) write(44,50) (tt*dstep(I),I=6,nn) cc TestStep to next stationary point call newton(hessian,grad,dstep,dnst,gmat,NN) write(44,*) ' dNewtonStep norm whould be: ', dnst c(bohr,radian)-->>(angstroem,degree) call bohrtoang(point,NN) call vec_dcopy(tang,told,nn) return 50 FORMAT(1X,3F20.10) end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC subroutine buildk(dk,req,tang,rgr,nn,pla) c build the K-Matrix real*8 dk,req,tang,rgr,pla integer nn,i,j dimension dk(nn,nn+1),req(nn-1,nn),tang(nn),rgr(nn-1) do 10,i=1,nn-1 do 11,j=1,nn dk(i,j)=req(i,j) 11 continue 10 continue do 20,j=1,nn dk(nn,j)=tang(j) 20 continue do 30,i=1,nn-1 dk(i,nn+1)=-rgr(i) 30 continue ccccc eventuell: ersetze pla=zero durch kleinen Schritt dk(nn,nn+1)=pla return end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine newton(hi,fi,step,stepl,gmat,nn) c calculates STEP = - HI^(-1) * FI and returns step and steplength real*8 hi,fi,step,stepl,gmat,gl,vec_mdot integer i,nn dimension hi(nn,nn),fi(nn),step(nn),gmat(nn,nn),gl(nn,nn+1) call mat_dcopy(hi,gl,nn,nn) do 10,i=1,nn gl(i,nn+1)=-fi(i) 10 continue call linsolve(gl,step,nn) stepl=dsqrt(vec_mdot(step,gmat,step,nn)) return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine mat_diag(d,n,v) c create a NxN diagonal matrix D with entries V real*8 d,v integer n,i dimension d(n,n) if(n .le. 0) return call vec_dinit(d,n*n,0.0d0) do i=1,n d(i,i)=v end do return end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine vec_dinit(arra1,ndim,wert) c Set all elements of ARRA1 to the double value of WERT real*8 arra1,wert integer ndim,i dimension arra1(ndim) if(ndim .le. 0) return do i=1,ndim arra1(i)=wert end do return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine vec_madd * (nr1,mat1,d1,fact1,nr2,mat2,d2,fact2,nrx,matx,dx,dim) c matx(nrx)=fact1*mat1(nr1)+fact2*mat2(nr2) c d1,d2,dx - columns c dim - rows real*8 mat1,mat2,matx,fact1,fact2,vec integer nr1,nr2,nrx,d1,d2,dx,dim,i dimension mat1(dim,d1),mat2(dim,d2),matx(dim,dx),vec(dim,1) call vec_dinit(vec,dim,0.0d0) do 10,i=1,dim vec(i,1)=fact1*mat1(i,nr1)+fact2*mat2(i,nr2) 10 continue call vec_mat_centry(vec,matx,nrx,dim,dx,'i') return end c ccccccccccccccccccccccccccccccccccccccccccccccccc subroutine vec_mat_centry(vec,mat,npos,ndim,ncol,zkind) c Copies VEC 'i'n or 'o'ut the NPOS column of matrix MAT real*8 vec,mat integer npos,ndim,ncol,i character zkind dimension vec(ndim),mat(ndim,ncol) if(npos.gt.ncol) return if (zkind.eq.'i') then do i=1,ndim mat(i,npos)=vec(i) end do endif if (zkind.eq.'o') then do i=1,ndim vec(i)=mat(i,npos) enddo endif return end c subroutine vec_mat_rentry(vec,dmat,npos,ndim,nrow,zkind) c Copies VEC 'i'n or 'o'ut the NPOS row of matrix MAT real*8 vec,dmat integer npos,ndim,nrow,i character zkind dimension vec(ndim),dmat(nrow,ndim) if(npos.gt.nrow) return if (zkind.eq.'i') then do 10,i=1,ndim 10 dmat(npos,i)=vec(i) endif if (zkind.eq.'o') then do 20,i=1,ndim 20 vec(i)=dmat(npos,i) endif return end c subroutine mat_dcopy(arra1,arra2,ndim,mdim) c Copies matrix ARRA1 to matrix ARRA2 real*8 arra1,arra2 integer i,j,ndim,mdim dimension arra1(ndim,mdim),arra2(ndim,mdim) do 10,i = 1,ndim do 20,j = 1,mdim arra2(i,j) = arra1(i,j) 20 continue 10 continue return end c real*8 function vec_mdot(vec1,metric,vec2,dim) c metric scalar product of VEC1 and VEC2 with respect to METRIC integer dim real*8 vec1,metric,vec2,cov2,vec_ddot dimension vec1(dim),vec2(dim),cov2(dim),metric(dim,dim) call matmult(metric,dim,vec2,1,cov2,dim) vec_mdot=vec_ddot(vec1,cov2,dim) return end c real*8 function vec_ddot(arra1,arra2,ndim) c cartesian scalar product of ARRA1 and ARRA2 real*8 arra1,arra2 integer ndim,i dimension arra1(ndim),arra2(ndim) vec_ddot = 0 if(ndim .le. 0)return do 70 i = 1,ndim vec_ddot = vec_ddot + arra1(i) * arra2(i) 70 continue return end c subroutine matmult(m1,d1,m2,d2,mout,dim) c matrix multiplication: m1[d1,dim]*m2[dim,d2]=mout[d1,d2] real*8 m1,m2,mout,swap integer d1,d2,dim,i,j,k dimension m1(d1,dim),m2(dim,d2),mout(d1,d2),swap(d1,d2) call vec_dinit(swap,d1*d2,0.0d0) do 10,i=1,d1 do 20,j=1,d2 do 30,k=1,dim swap(i,j)=swap(i,j)+m1(i,k)*m2(k,j) 30 continue 20 continue 10 continue call vec_dcopy(swap,mout,d1*d2) return end c subroutine linsolve(matein,x,n) c solves MATEIN * X = 0 real*8 matein,dre,x,w integer n dimension matein(n,n+1),dre(n,n+1),x(n),w(n) call householder(matein,dre,w,n,n+1) call loesung(dre,x,n) return end c subroutine householder(matein,mataus,w,n,sp) c householder transformation integer k,l,m,n,sp real*8 alpha,rho,sk,s,w,matein,mataus dimension matein(n,sp),mataus(n,sp),w(n) call vec_dcopy(matein,mataus,n*sp) do 10,k=1,n-1 s=0.0d0 do 20,l=k,n s=s+mataus(l,k)*mataus(l,k) 20 continue if (mataus(k,k).lt.0) then alpha=dsqrt(s) else alpha=-dsqrt(s) endif rho=dsqrt(s-mataus(k,k)*mataus(k,k)+ * (mataus(k,k)-alpha)*(mataus(k,k)-alpha)) if(rho .lt. 1.0d-63) then rho= 1.0d-63 C w(k)= 0.0d0 C goto 22 endif w(k)=(mataus(k,k)-alpha)/rho 22 continue do 30,l=k+1,n w(l)=mataus(l,k)/rho 30 continue do 40,l=k,sp sk=0.0d0 do 55,m=k,n sk=sk+w(m)*mataus(m,l) 55 continue do 60,m=k,n mataus(m,l)=mataus(m,l)-2*sk*w(m) 60 continue 40 continue 10 continue return end c subroutine loesung(matein,x,n) c intern linear equation routine real*8 matein,x integer i,j,n,k dimension matein(n,n+1),x(n) do 10,i=1,n 10 x(i)=matein(i,n+1) x(n)=x(n)/matein(n,n) do 20,i=1,n-1 k=n-i do 30,j=k+1,n 30 x(k)=x(k)-matein(k,j)*x(j) x(k)=x(k)/matein(k,k) 20 continue return end c subroutine dfp(einmat,q,p,nn,upmat) ccc DFP-Update of Hessian ccc Heidrich,Kliesch,Quapp: Properties... ccc S.51 Eqn.(21)/(22) ccc ccc einmat Hessian to update (H_0) ccc q difference of gradient (g_1 - g_0) ccc p difference of point (x_1 - x_0) ccc upmat updated Hessian (H_1) real*8 einmat,q,p,upmat,qp,unit,b1,b2,b3, 1 vec_ddot,eps,sw integer nn dimension unit(nn,nn),b3(nn,nn), + einmat(nn,nn),upmat(nn,nn), + b1(nn,nn),b2(nn,nn),sw(nn,nn), + q(nn), p(nn) data eps/1.0d-30/ call vec_dinit(upmat,nn*nn,0.0d0) call mat_diag(unit,nn,1.0d0) qp=vec_ddot(q,p,nn) ccc if(dabs(qp).lt.eps) call cerr(-31) call matmult(q,nn,p,nn,sw,1) call matadd(unit,(-1.0d0)/qp,sw,b1,nn,nn) call matmult(p,nn,q,nn,sw,1) call matadd(unit,(-1.0d0)/qp,sw,b2,nn,nn) call matmult(q,nn,q,nn,b3,1) call tri_mat(b1,nn,nn,einmat,nn,b2,nn,sw) call matadd(sw,1.0d0/qp,b3,upmat,nn,nn) return end c subroutine tri_mat(dm1,nr1,nc1,dm2,nc2,dm3,nc3,dres) c dm1[nr1,nc1]*dm2[nc1,nc2]*dm3[nc2,nc3]=dres[nr1,nc3] real*8 dm1,dm2,dm3,dres,r1 integer nr1,nc1,nc2,nc3 dimension dm1(nr1,nc1),dm2(nc1,nc2),dm3(nc2,nc3), + r1(nr1,nc2),dres(nr1,nc3) call matmult(dm1,nr1,dm2,nc2,r1,nc1) call matmult(r1,nr1,dm3,nc3,dres,nc2) return end c subroutine ortcompl(vec,compl,dmet,ndim) c calculate an Orthonormalcomplement COMPL to vector VEC real*8 basis,vec,compl,swap,det,dmgs,dmet integer i,j,ndim dimension vec(ndim),compl(ndim-1,ndim),basis(ndim,ndim), * swap(ndim),dmet(ndim,ndim) call vec_dinit(basis,ndim*ndim,0.0d0) do 10,i=1,ndim basis(i,i)=1.0d0 10 continue do 20,i=1,ndim if (vec(i).ne.0.0d0) then do 30,j=1,ndim swap(j)=basis(j,1) basis(j,1)=vec(j) if (i.ne.1) basis(j,i)=swap(j) 30 continue goto 40 endif 20 continue 40 continue det=dmgs(basis,dmet,ndim) call vec_dinit(compl,(ndim-1)*ndim,0.0d0) call mat_trans(basis,ndim,ndim,basis) do 60,i=2,ndim call vec_mrowcopy(i,basis,ndim,i-1,compl,ndim-1,ndim) 60 continue return end c real*8 function dmgs(a,dmet,dim) c modified Gram-Schmidt algorithm return det(a) real*8 det,rr,a,q,r,vec_mmdot,dmet,dzero integer dim,i,k dimension a(dim,dim),q(dim,dim),r(dim,dim),dmet(dim,dim) data dzero/0.0d0/ call vec_dinit(q,dim*dim,0.0d0) call vec_dinit(r,dim*dim,0.0d0) det=1.0d0 call mat_dcopy(a,q,dim,dim) do 10,k=1,dim if (k.eq.1) goto 21 do 20,i=1,k-1 r(i,k)=vec_mmdot(i,q,dim,k,q,dim,dim,dmet) call vec_madd(k,q,dim,1.0d0,i,q,dim,-r(i,k),k,q,dim,dim) 20 continue 21 continue rr=vec_mmdot(k,q,dim,k,q,dim,dim,dmet) r(k,k)=dsqrt(rr) det=det*r(k,k) if(rr.lt.dzero) goto 33 call vec_madd(k,q,dim,1.0d0/r(k,k),k,q,dim,0.0d0,k,q,dim,dim) 10 continue call mat_dcopy(q,a,dim,dim) dmgs=det return 33 CONTINUE cc call cerr(-23) return end c subroutine bof(einmat,q,p,nn,upmat) cc Bofill Update of Hessian, see cc JCC 15 (1994) 1-11, eqs.(10)-(13) cccccccccccccccccccccccccccccccccccccccccccccccccccc cc einmat Hessian to update (H_0) c cc q (=gamma) difference of gradient (g_1 - g_0) c cc p (=delta) difference of point (x_1 - x_0) c cc upmat updated Hessian (H_1) c cc zi q - (H_0) p c cc Hms Murtag-Sargent update ccc cc Hpo Powell update c cc phi optimal combination between Hms/Hpo c cc c cccccccccccccccccccccccccccccccccccccccccccccccccccccc real*8 einmat,q,p,upmat,zi,unit,b1,v2, 1 vec_ddot,eps,sw,pzi,phi,pp,zizi,Hms,Hpo integer nn dimension unit(nn,nn), + einmat(nn,nn),upmat(nn,nn), + b1(nn,nn),sw(nn,nn),v2(nn), + q(nn),p(nn),zi(nn),Hms(nn,nn),Hpo(nn,nn) data eps/1.0d-30/ call vec_dinit(upmat,nn*nn,0.0d0) call vec_dinit(Hms,nn*nn,0.0d0) call vec_dinit(Hpo,nn*nn,0.0d0) call mat_diag(unit,nn,1.0d0) call matmult(einmat,nn,p,1,v2,nn) call vec_dmult_add(q,v2,nn,-1.0d0,zi) Cc Hms Murtag-Sargent update -- eq.(11) pzi=vec_ddot(p,zi,nn) ccc if(dabs(pzi).lt.eps) call cerr(-31) call matmult(zi,nn,zi,nn,b1,1) call matadd(einmat,1.0d0/pzi,b1,Hms,nn,nn) Cc Hpo Powell update -- eq.(12) pp=vec_ddot(p,p,nn) cccc if(dabs(pp).lt.eps) call cerr(-31) zizi=vec_ddot(zi,zi,nn) cccc if(dabs(zizi).lt.eps) call cerr(-31) call matmult(p,nn,p,nn,sw,1) call matadd(einmat,(-1.0d0)*pzi*pzi/pp/pp,sw,Hpo,nn,nn) call matmult(zi,nn,p,nn,b1,1) call matadd(Hpo,1.0d0/pp,b1,sw,nn,nn) call matmult(p,nn,zi,nn,b1,1) call matadd(sw,1.0d0/pp,b1,Hpo,nn,nn) c c combination of Hms and Hpo by phi -- eqs.(10,13) phi= 1.0d0 - pzi*pzi/pp/zizi call vec_dinit(sw,nn*nn,0.0d0) call matadd(sw,(1.0d0-phi),Hms,b1,nn,nn) call matadd(b1,phi,Hpo,upmat,nn,nn) 2 continue return end subroutine gettang(curveq,gmat,tangent,nn) c This routine calculate the tangent on a curve c whereas the tangent is unique determined by: c c H*t=0 c |t|=1 c det[H,t]>0 c real*8 tangent,gmat,curveq,q,ss1,det,zero,dmnorm integer nn,i,j dimension tangent(nn),curveq(nn-1,nn),q(nn*nn),ss1(nn,nn) c dimension tangent(nn),curveq(nn-1,nn),q(nn,nn),ss1(nn,nn) dimension gmat(nn,nn) data zero /0.0d0/ call qr(curveq,q,nn-1) call vec_dcopy(q(nn*(nn-1)+1),tangent(1),nn) c do 1,i=1,nn c 1 tangent(i)=q(nn,i) c wq call vec_dmult_constant + (tangent,nn,1/dmnorm(tangent,gmat,nn),tangent) call mat_dcopy(curveq,ss1,nn-1,nn) do 10,i=1,nn 10 ss1(nn,i)=tangent(i) if (det(ss1,nn).lt.zero) then call vec_dmult_constant(tangent,nn,-1.0d0,tangent) endif return end c subroutine qr(gex,q,n) C QR-Zerlegung der Jacobimatrix von ge (gex) integer i,j,k,n real*8 s1,s2,s,gex,q,r,qq,rr dimension gex(n,n+1),q(n+1,n+1),r(n+1,n), * qq(n+1,n+1),rr(n+1,n) do 10,i=1,n do 20,j=1,n+1 r(j,i)=gex(i,j) 20 continue 10 continue do 30,i=1,n+1 do 40,j=1,n+1 q(i,j)=0.0d0 40 continue q(i,i)=1.0d0 30 continue do 50,i=1,n do 60,j=i+1,n+1 s1=r(i,i) s2=r(j,i) s=dsqrt(s1*s1+s2*s2) if (s.gt.0.0d0) then s1=s1/s s2=s2/s else s1=0.0d0 s2=0.0d0 endif call vec_dcopy(r,rr,(n+1)*n) do 70,k=1,n rr(i,k)=s1*r(i,k)+s2*r(j,k) rr(j,k)=s1*r(j,k)-s2*r(i,k) 70 continue call vec_dcopy(rr,r,(n+1)*n) call vec_dcopy(q,qq,(n+1)*(n+1)) do 80,k=1,n+1 qq(i,k)=s1*q(i,k)+s2*q(j,k) qq(j,k)=s1*q(j,k)-s2*q(i,k) 80 continue call vec_dcopy(qq,q,(n+1)*(n+1)) 60 continue 50 continue do 90,i=1,n+1 do 100,j=1,n+1 q(i,j)=qq(j,i) 100 continue 90 continue return end c c subroutine recorrect(corr,cold,gmat,ocorr,nn) c manipulation of the corrector for a better convergence c integer nn c real*8 corr,cold,dmnorm,gmat,vec_mdot c logical olen,odir,ocorr,odone c dimension corr(nn),cold(nn),gmat(nn,nn) c ocorr=.false. c if (dmnorm(corr,gmat,nn).gt.dmnorm(cold,gmat,nn)) then c olen=.true. c else c olen=.false. c endif c if (vec_mdot(corr,gmat,cold,nn).lt.0.0d0) then c odir=.true. c else c odir=.false. c endif c odone=.false. c if (odir.and.olen.and.ocorr) then c call vec_dmult_constant(corr,nn,0.6d0,corr) c odone=.true. c endif c ocorr=odone c return c end c real*8 function dmnorm(vec,gm,nn) c metric norm of VEC with resp. to metric GM real*8 vec,vec_mdot,gm integer nn dimension vec(nn),gm(nn,nn) dmnorm=dsqrt(vec_mdot(vec,gm,vec,nn)) return end c real*8 function det(matein,nn) c calculation of the determinant of MATEIN integer i,nn real*8 matein,mataus,w dimension matein(nn,nn),mataus(nn,nn),w(nn) call householder(matein,mataus,w,nn,nn) det=1.0d0 do 10,i=1,nn if (dabs(mataus(i,i)).lt.1.0d-25) then det=0.0d0 goto 20 endif if (dabs(det) .lt.1.0d-25) then det=0.0d0 goto 20 endif 10 det=det*mataus(i,i) if (iand(nn,1).eq.0) * det=-det 20 continue return end c real*8 function adet(mat,sdet,dim) c calculate the Determinant of MAT real*8 mat,amat,det,fact,sdet,sz integer dim,itest,sp,nexp,k,dreieck,ll dimension mat(dim,dim),amat(dim,dim),sp(dim-1),fact(dim) sz=1.0d0/16.0d0 call mat_dcopy(mat,amat,dim,dim) itest=dreieck(amat,sp,dim) do 5,k=1,dim-1 5 continue if (itest.eq.-1) then adet=0.0d0 sdet=0.0d0 else nexp=0 det=1 sdet=1.0d0 do 10,k=1,dim det=det*amat(k,k) fact(k)=amat(k,k) if (k.lt.dim) then if (k.lt.sp(k)) det=-det endif 20 continue if (dabs(det).lt.sz) then det=det*16.0d0 nexp=nexp-1 goto 20 endif 30 continue if (dabs(det).ge.1.0d0) then det=det*sz nexp=nexp+1 goto 30 endif if (k.eq.dim-1) then sdet=det ll=nexp endif 10 continue if (nexp.lt.0) then nexp=-nexp do 40,k=1,nexp det=det*sz 40 continue elseif (nexp.gt.0) then do 50,k=1,nexp det=det*16.0d0 50 continue endif if (ll.lt.0) then ll=-ll do 410,k=1,ll sdet=sdet*sz 410 continue elseif (ll.gt.0) then do 510,k=1,ll sdet=sdet*16.0d0 510 continue endif adet=det endif return end c subroutine mat_trans(matin,nrow,ncol,matout) c transpose MATOUT of matrix MATIN real*8 matin,matout,swap integer nrow,ncol,i,j dimension matin(nrow,ncol),matout(ncol,nrow),swap(ncol,nrow) do 10,i=1,nrow do 20,j=1,ncol swap(j,i)=matin(i,j) 20 continue 10 continue call mat_dcopy(swap,matout,ncol,nrow) return end cccccccccccccccccccccccccccccccccccccccccccccc c subroutine toang(vec1,vec2,nn) cc (bohr,radian)->(angstroem,degree) for internal coordinates c real*8 vec1,vec2,torad,btoa c integer i,nr,nn c dimension vec1(nn),vec2(nn) c data btoa /0.52917706d0/ c torad=dacos(-1.0d0)/180.0d0 c nr=(nn+6)/3-1 c do i=1,nn c if (i.le.nr) then c vec2(i)=vec1(i)*btoa c else c vec2(i)=vec1(i)/torad c endif c end do c return c end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine angtobohr(point,NN) real*8 torad,btoa,point(nn) integer i,nn,n3 c (bohr,radian) <<-- (angstroem,degree) for internal coordinates c torad=dacos(-1.0d0)/180.0d0 torad=3.141592653589793d0/180.0d0 btoa =0.52917706d0 n3=(nn+6 +2 )/3 c c point(1)=point(1)/btoa c point(2)=point(2)/btoa c point(3)=point(3)*torad c If(nn.gt.3) then c Do 1 I=1,N3-3 c point(3*I+1)=point(3*I+1)/btoa c point(3*I+2)=point(3*I+2)*torad c1 point(3*I+3)=point(3*I+3)*torad c Endif c note: fix dih4 , reorganize the list of coordinates point(1)=point(1)/btoa point(2)=point(2)/btoa point(3)=point(3)*torad point(4)=point(4)/btoa point(5)=point(5)*torad If(nn.gt.5) then Do 1 I=2,N3-4 point(3*I)=point(3*I)/btoa point(3*I+1)=point(3*I+1)*torad 1 point(3*I+2)=point(3*I+2)*torad Endif point(nn-1)=point(nn-1)/btoa point(nn)=point(nn)*torad return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine bohrtoang(point,NN) real*8 torad,btoa,point(nn) integer i,nn,n3 c (bohr,radian)-->>(angstroem,degree) for internal coordinates c torad=dacos(-1.0d0)/180.0d0 torad=3.141592653589793d0/180.0d0 btoa =0.52917706d0 n3=(nn+6 +2 )/3 c note dih6 is fixed, also dih22 point(1)=point(1)*btoa point(2)=point(2)*btoa point(3)=point(3)/torad point(4)=point(4)*btoa point(5)=point(5)/torad If(nn.gt.5) then Do 1 I=2,N3-4 point(3*I)=point(3*I)*btoa point(3*I+1)=point(3*I+1)/torad 1 point(3*I+2)=point(3*I+2)/torad Endif point(nn-1)=point(nn-1)*btoa point(nn)=point(nn)/torad return end ccccccccccccccccccccccccccccccccccccccccccccccccc subroutine vec_dmult_add(arra1,arra2,ndim,zahl,arra4) c For each Komponent: ARRA4 = ARRA1 + ZAHL*ARRA2 real*8 arra1,arra2,arra4,zahl integer ndim,i dimension arra1(ndim),arra2(ndim),arra4(ndim) if(ndim .le. 0)return do 145 i = 1,ndim arra4(i) = arra1(i)+arra2(i)*zahl 145 continue return end cccccccccccccccccccccccccccccccccccccccccccccc subroutine vec_dcopy(arra1,arra2,ndim) c Copies ARRA1 to ARRA2 real*8 arra1,arra2 integer i,ndim dimension arra1(ndim),arra2(ndim) if(ndim .le. 0) return do 10 i = 1,ndim arra2(i) = arra1(i) 10 continue return end ccccccccccccccccccccccccccccccccccccccccccccccccc subroutine vec_dmult_constant(arra1,ndim,zahl,arra2) c ARRA2 = ZAHL * ARRA1 real*8 arra1,arra2,zahl integer ndim,i dimension arra1(ndim),arra2(ndim) if(ndim .le. 0)return do i = 1,ndim arra2(i) = arra1(i)*zahl end do return end ccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine vec_mrowcopy(nri,matin,din,nra,matout,dout,dim) c copies the (nri) row of (matin) with (din) row and (dim) col c to the (nra) row of (matout) with (dout) row and (dim) col real*8 matin,matout integer nri,din,nra,dout,dim,i dimension matin(din,dim),matout(dout,dim) ccc if ((nri.gt.din).or.(nra.gt.dout)) call cerr(17) do 10,i=1,dim matout(nra,i)=matin(nri,i) 10 continue return end cccccccccccccccccccccccccccccccccccccccccccccccccc real*8 function vec_mmdot &(nri,matin,din,nra,matout,dout,dim,dmet) c multiply the (nri) col of (matin) and the (nra) col of c (matout). (dim) is the number of the rows with metric dmet real*8 matin,matout,prod,v1,v2,dmet,vec_mdot integer nri,din,nra,dout,dim dimension matin(dim,din),matout(dim,dout),dmet(dim,dim), + v1(dim),v2(dim) ccc if ((nri.gt.din).or.(nra.gt.dout)) call cerr(-100) call vec_mat_centry(v1,matin,nri,din,dim,'o') call vec_mat_centry(v2,matout,nra,dout,dim,'o') prod=vec_mdot(v1,dmet,v2,dim) vec_mmdot=prod return end ccccccccccccccccccccccccccccccccccccccccccccccc subroutine matadd(dm1,fact,dm2,out,n,m) c OUT(N,M) = DM1(N,M) + FACT * DM1(N,M) real*8 dm1,fact,dm2,out integer m,n dimension dm1(n,m),dm2(n,m),out(n,m) call vec_dmult_add(dm1,dm2,n*m,fact,out) return end cccccccccccccccccccccccccccccccccccccccccccc integer function dreieck(mat,sp,dim) c Dreiecksfaktorisierung of MAT ; dim(SP) = DIM-1 real*8 mat,z integer sp,dim,it,k,i,s,j dimension mat(dim,dim),sp(dim-1) it=0 do 10,k=1,dim-1 z=0.0d0 do 20,i=k,dim if (dabs(mat(i,k)).gt.z) then z=dabs(mat(i,k)) s=i endif 20 continue if (z.eq.0.0d0) then dreieck=-1 return endif sp(k)=s if (k.lt.s) then do 30,j=1,dim z=mat(k,j) mat(k,j)=mat(s,j) mat(s,j)=z 30 continue endif do 40,i=k+1,dim mat(i,k)=mat(i,k)/mat(k,k) 40 continue do 50,j=k+1,dim do 60,i=k+1,dim mat(i,j)=mat(i,j)-mat(i,k)*mat(k,j) 60 continue 50 continue 10 continue if (mat(dim,dim).eq.0.0d0) it=-1 dreieck=it return end cccccccccccccccccccccccc