PROGRAM HCNvri2 ! ! par t(2) rgf to dir=rsel ! ! W.Quapp, Benjamin Schmidt, Decembre 2009 ! ! nn = Dimension, N3 = number of atoms ! le = Chainlength (even integer) ! iterNTmax= maximal number of NT correctors: if corector ! uses more than, say 4, steps: it works not fine.. ! pst,pstl = Steplength of predictor, depends from iteration ! eps = Convergence in RGF corrector, red from part 1 ! is useless for predictors only ! kette = Chain between points xstart and xend, ! xend is useless for predictors only ! rsel = Searchdirection RGF method = rsel in RGF corrector ! ! scriHCNntVRI: is the working code in Linux, combining the ! followings FORTRAN programms, which use following steps ! (1) Start: Initial settings, with ! (i) Fix start- and end point of the chain ! (ii) Fix searchdirection "dir", and initial guess of "VRI" ! (iii) Put straigt chain of "length" between both points ! ! this part: ! (2) To "dir" calculate corresponding NT, by pred-corr, ! or by predictors along AG only ! ! lather parts: ! (3) Calculate VRI approximation, ! one version: use its gradient for new "dir" ! other version: use only Ag for a predictor step ! (4) Organize points for Gamess-US input for (2) or (3) ! HCNvriReadOUT: translates output into usable input files ! (5) Organize Iterations up to convergence of VRI ! by repeating of (2)-(4) ! ! ! INTEGER nn,N3,le,iterNTmax DOUBLE PRECISION stl,pstl PARAMETER (nn=3,N3=3,le=50,iterNTmax=10,stl=1.d0,pst=0.012d0) ! ! RGF-Verfahren ! DOUBLE PRECISION kette(le+1,nn),rkette(le+1,nn),eps,tgt,tt,pla, 1 rsel(nn),rc(nn),grad(nn),H(nn,nn),Proj(nn,nn),Ag(nn),tmin, 2 xstart(nn),xend(nn),sw(nn),prestep(nn),Adj(nn,nn),Agmin, 3 norm,gmat(nn,nn),dstep(nn),detk,tminNT,AgmiNT,factor, 4 projNm1(nn-1,nn),dk(nn,nn+1),rgeq(nn-1,nn),Help(nn,nn), 5 rgrad(nn-1),tcov(nn),tang(nn),told(nn), 6 unit(nn,nn),dmnorm,vec_mdot,vec_ddot,det, 7 t1old(nn),t1(nn),rgnorm,point(nn),eival(nn),eivec(nn,nn), 8 gradcont(nn),ginv(nn,nn),energy,sgn INTEGER k,j,i,jnr,istatus,icalc,iter,iterNT,Natoms(N3) logical omat_inv Character*2 CharAtoms(N3),Zcoor(nn) Character*70 zmatLines(20) Data Zcoor/'hc','nc','w '/ Data CharAtoms/'H ','C ','N '/ Data Natoms/1,6,7/ c pure and applied chemistry, 51, 1 (1979): c angstroms per bohr data toang /0.52917706d0/ ccccccccccccccccccccccccccccccccccccccccc OPEN(8,FILE='hcnVkette.dat') c OPEN(91,FILE='hcnVkette.mol') c OPEN(95,FILE='hcnVkette.gau') ccccccccccccccccccccccccccccccccccccccccc OPEN(7,FILE='chain.dat',status='unknown',access='append') OPEN(9,FILE='hcnVketNrj.dat') OPEN(11,FILE='hcnVdir.dat') OPEN(13,FILE='hcnVener.dat') OPEN(53,FILE='hcnVEall.dat',status='unknown',access='append') OPEN(14,FILE='hcnVgrad.dat') OPEN(15,FILE='hcnVhesse.dat') OPEN(16,FILE='hcnVeps.dat') OPEN(17,FILE='hcnViter.dat') OPEN(18,FILE='hcnVistat.dat') OPEN(19,FILE='hcnVicalc.dat') OPEN(20,FILE='hcnVketst.dat') OPEN(21,FILE='hcnVketen.dat') OPEN(25,FILE='hcnVtold.dat') OPEN(26,FILE='hcnVt1old.dat') OPEN(28,FILE='hcnVketNrj1.dat') OPEN(32,FILE='hcnVnode.dat') OPEN(33,FILE='hcnVmin.dat') OPEN(48,FILE='hcnVgmatr.dat') c protocol file: all interesting output of the vri-search OPEN(44,FILE='hcnVprotocol.txt',status='unknown',access='append') cccccccccccccccccccccccccccccccccccccccc 2 FORMAT(3F20.14) N=nn c j1=le/2 j1=1 k1=1 c k1=le+1 c j1,k1 for iteration in part 3 rewind 28 write(28,*) j1,k1 rewind 16 read(16,*) eps c eps is criterion for rgf convergence to given rsel-NT rewind 9 Read(9,*) jnr c jnr is the currend node number If(jnr.eq.2) then write(*,*) 'eps = ',eps write(44,*)'eps = ',eps endif rewind 8 DO 1 k=1,le+1 Read(8,*) (kette(k,i),i=1,nn) 1 CONTINUE rewind 19 read(19,*) icalc,sgn rewind 17 Read(17,*) iter,iterNT if(iter.eq.0) then pstl= pst else pstl= 2.0d0*pst endif if(iterNT.gt.iterNTmax) then write(*,*) write(*,*)'Stop: NT-node not convergent at jnr',jnr 1 ,' in Iteration',iter write(*,*) write(44,*) write(44,*)'Stop: NT-node not convergent at jnr',jnr 1 ,' in Iteration',iter write(44,*) if(jnr.lt.le+1) then jnr=jnr+1 else istatus=0 rewind 18 write(18,*) istatus sgn=-1.0d0 icalc=0 rewind 19 write(19,*) icalc,sgn goto 66 endif iterNT=0 write(17,*) iter,iterNT endif if(icalc.eq.1) then write(44,*) 'kette(',jnr,') in angstroem,degree HCNvri2 ' Do 7, i=1,nn 7 write(44,*) Zcoor(i), kette(jnr,i) endif c c ccccccccccccccccccccccccccccccccccccccccccccccccc c (bohr,radian)<->(angstroem,degree) for internal coordinates c current point (angstroem,degree) ->(bohr,radian) torad=dacos(-1.0d0)/180.0d0 btoa =0.52917706d0 point(1)=kette(jnr,1)/btoa point(2)=kette(jnr,2)/btoa point(3)=kette(jnr,3)*torad c Note: c in corrector loop with use of gradient and hessian, there c are used bohr and radian. Input of Gamess, however, is in c angstroem and degree !!! c c ccccccccccccccccccccccccccccccccccccccccccccccccc call mat_diag(unit,nn,1.0d0) rewind 11 Read(11,*) (rsel(i),i=1,nn) if(jnr.eq.2) then write(*,*) 'dir in vri2 ', (rsel(i),i=1,nn) write(44,*) 'dir in vri2 ', (rsel(i),i=1,nn) endif rewind 18 read(18,*) istatus c if(istatus.eq.10) then write(44,*)'Abbruch in vri2 cycle' write(*,*) 'Abbruch in vri2 cycle' goto 66 else istatus=2 rewind 18 write(18,*) istatus endif c rewind 25 Read(25,*) (told(i),i=1,nn) rewind 26 Read(26,*) (t1old(i),i=1,nn) c rewind 19 read(19,*) icalc,sgn if(icalc.eq.0) then icalc=1 write(19,*) icalc goto 66 else icalc=0 write(19,*) icalc rewind 13 Read(13,*) energy write(44,*)'Energy ',energy,' at iter ',iter, 1 'kete(',jnr,') iterNT=',iterNT write(53,*)'Energy ',energy,' at iter ',iter, 1 'kete(',jnr,') iterNT=',iterNT rewind 14 Read(14,*) (grad(i),i=1,nn) rewind 15 Do 15 k=1,nn 15 Read(15,*) (H(k,i),i=1,nn) endif cc use 'Proj' for a help matrix c call eigen(H,nn,nn,eival,eivec,Proj,0) c write(44,*)'eigenvalues of H ' c write(44,*)(eival(i),i=1,nn) c write(44,*)'eigenvectors of H ' c do k=1,nn c write(44,*)(eivec(k,i),i=1,nn) c enddo c c ccccccccccccccccccccccccccccccccccccccccccccccccc c Ginv is the G-Matrix of Wilson: Ginv=B*B^T=g^{ij} c Gmat=Ginv^{-1} is the usual covarinat metric matrix g_{ij} c note: all is without mass-weighting rewind 48 Do 766, k=1,nn 766 READ(48,*) (ginv(k,I),I=1,nn) READ(48,*) Do 767, k=1,nn 767 READ(48,*) (gmat(k,I),I=1,nn) READ(48,*) IF(jnr.eq.1) then write(44,*) ' ginv is from the B-mat ' Do 7711, k=1,nn 7711 WRITE(44,2) (ginv(k,I),I=1,nn) write(44,*) ' gmat ' Do 7712, k=1,nn 7712 WRITE(44,2) (gmat(k,I),I=1,nn) Endif c testprint c CALL matmult(ginv,nn,gmat,nn,unit,nn) c write(*,*)'test metric ginv*gmat' c do k=1,nn c write(*,*)(unit(k,i),i=1,nn) c enddo c unit, gmat, ginv is for Metric in internal coordinates call vec_dinit(gradcont,nn,0.0d0) call matmult(ginv,nn,grad,1,gradcont,nn) write(44,*) ' grad ' WRITE(44,2) (grad(I),I=1,nn) write(44,*) ' gradcont ' WRITE(44,2) (gradcont(I),I=1,nn) gnorm=dsqrt(vec_ddot(grad,gradcont,nn)) WRITE(44,*) ' PES full gradient norm ',gnorm ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c rsel has to be covariant, like the gradient c projNm1: n-1 linearly independent lines from Proj call matmult(ginv,nn,rsel,1,rc,nn) c write(44,*) 'dir contravariant ',(rc(i),i=1,nn) CALL projector(Proj,rsel,rc,nn) c old, is ruled out: call ortcompl(rsel,projNm1,unit,nn) c write(44,*) ' Proj ' c Do 7713, k=1,nn c7713 WRITE(44,2) (Proj(k,I),I=1,nn) do 7719 k=1,nn-1 do i=1,nn projNm1(k,i)=Proj(k,i) enddo 7719 continue cccccccccccccccccccccccccccccccccccccccccccccccccccc c test print : is o.k. c write(44,*) ' rsel ' c WRITE(44,*) (rsel(I),I=1,nn) c call matmult(proj,nn,rsel,1,sw,nn) c write(44,*) ' projNm1 * rsel ' c WRITE(44,*) (sw(I),I=1,nn) cccccccccc WQ 2009 cccccccccccccccccccccccccccccccccccccc c c c proj: 1 x covariant, 1 x contravariant tensor c c 'Reducting' of Hessian: is again 2x covariant c c rgeq is again two-fold covariant c c c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc call matmult(projNm1,nn-1,H,nn,rgeq,nn) write(44,*) ' rgeq ' Do 7714, k=1,nn-1 7714 WRITE(44,2) (rgeq(k,I),I=1,nn) call vec_dinit(rgrad,nn-1,0.0d0) call matmult(projNm1,nn-1,grad,1,rgrad,nn) c sw: use full dimension of rgrad for metric call matmult(Proj,nn,grad,1,sw,nn) write(44,*) ' rgrad ' WRITE(44,*) (sw(I),I=1,nn) rgnorm=dmnorm(sw,ginv,nn) WRITE(44,*) ' PES reduced grad norm ',rgnorm c Tangent ccccccccccccccccccccccccccccccccccccccccccccc call gettang(rgeq,gmat,tang,nn) write(44,*) 'tang from rgeq used ' WRITE(44,*) (tang(I),I=1,nn) c near VRI it is better to use Ag for the tangent c CALL mat_adj(H,nn,Adj) c CALL matmult(Adj,nn,grad,1,sw,nn) c tgt=dmnorm(sw,gmat,nn) c do i=1,nn c sw(i)=sw(i)/(-1.0d0*tgt) c enddo c write(44,*) 'tang by Ag' c WRITE(44,*) (sw(I),I=1,nn) ccccccccccccccccccccccccccccccccccccccccccccccccccccccc call matmult(rgeq,nn-1,tang,1,sw,nn) write(44,*) 'rgeq*tang gives ',(sw(I),I=1,nn-1) tgt=vec_mdot(tang,gmat,told,nn) write(44,*)'tang*G_ij*told ', tgt call vec_dinit(tcov,nn,0.0d0) call matmult(gmat,nn,tang,1,tcov,nn) write(44,*) ' tang is ' write(44,*) (tang(I),I=1,nn) write(44,*) ' tcov is ' write(44,2) (tcov(I),I=1,nn) cccccccccccccccccccccccccccccccccccccccccccccc call buildk(dk,rgeq,tcov,rgrad,nn) cccccccccccccccccccccccccccccccccccccccccccccc write(44,*) ' K mat is ' Do 7715, k=1,nn 7715 WRITE(44,*) (dk(k,I),I=1,nn+1) detk=det(dk,nn) write(*,*) ' detK is ', detk write(44,*) ' detK is ', detk if (detk.lt.0.0d0) then call vec_dmult_constant(tang,nn,-1.0d0,t1) else call vec_dcopy(tang,t1,nn) endif tgt=vec_mdot(tang,gmat,told,nn) if (tgt.lt.0.0d0) then call vec_dmult_constant(tang,nn,-1.0d0,tang) call vec_dmult_constant(tcov,nn,-1.0d0,tcov) call buildk(dk,rgeq,tcov,rgrad,nn) endif cccccccccccccccccccccccccccccccccccccccccccccc tg=vec_ddot(t1,t1old,nn) if (tg.lt.0.0d0 .and. jnr.gt.1 ) then write(*,*) ' bif of branches' write(44,*) ' bif of branch c that is ruled out ?? c call vec_dmult_constant(tang,nn,-1.0d0,tang) c call vec_dmult_constant(tcov,nn,-1.0d0,tcov) c call buildk(dk,rgeq,tcov,rgrad,nn) endif call linsolve(dk,dstep,nn) WRITE(44,*) 'dstep by linsolve would be ' WRITE(44,*)(dstep(i),i=1,nn) c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c ccccccc 1.order goto: 333 continue ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccccccccccccccccc Abbruch c Abbruch -- condition start c IF(jnr.gt.4 ) then write(44,*) ' Predictor steps only ************* pstl=',pstl write(*,*) ' Predictor steps only ************* pstl=',pstl c endif c IF( cc 1 Dabs(detk).lt.0.0005d0 .or. c 1 rgnorm.lt.eps .or. c 2 jnr.gt.4 .or. c 3 iterNT.gt.10 ) THEN ccc write(44,*) ' Predictor steps only **************** ' ccc write(*,*) ' Predictor steps only **************** ' write(44,*) 'Nr j in vri2 is =',jnr write(*,*) ' |rgnorm| ', rgnorm write(44,*) ' |rgnorm| ', rgnorm iterNT=1 rewind 17 write(17,*) iter,iterNT rewind 25 write(25,2) (tang(i),i=1,nn) rewind 26 write(26,2) (t1(i),i=1,nn) ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc if(jnr.lt.le+1) THEN cc do predictorstep point(1)=kette(jnr,1)/btoa point(2)=kette(jnr,2)/btoa point(3)=kette(jnr,3)*torad cc test predictor step by RGF 3 dk(nn,nn+1)=pstl call linsolve(dk,dstep,nn) tgt=dmnorm(dstep,gmat,nn) write(*,*) write(44,*)'pred-step length',tgt write(44,*)'pred-step by RGF3 would be' write(44,*)(dstep(i),i=1,nn) C C C C c C C C C c C C C C c C C C C c predictor step by Branin differential equ rewind 33 read(33,*) tmin, Agmin, AgmiNT CALL mat_adj(H,nn,Adj) CALL matmult(Adj,nn,grad,1,Ag,nn) tminNT=dmnorm(Ag,gmat,nn) write(44,*) 'tminNT of Ag ',tminNT factor=pstl/tminNT tgt=vec_mdot(Ag,gmat,told,nn) IF(jnr.lt.30 ) then if(tgt.lt. 0.0d0) factor=-1.0d0*factor endif IF(jnr.ge.30 )factor=-1.0d0*factor call vec_dmult_constant(Ag,nn,factor,dstep) tgt=dmnorm(dstep,gmat,nn) write(44,*)'pred-step length',tgt write(44,*)'pred-step by Branin is' write(44,*)(dstep(i),i=1,nn) c additional: search lowest node for |Ag| IF(tminNT .LT. AgmiNT ) THEN rewind 33 write(33,*) tmin, Agmin, tminNT ENDIF write(*,*)'along NT |Ag| = ',tminNT 1 ,' at node ', jnr write(44,*)'along NT |Ag| = ',tminNT write(44,*)' at node ', jnr c C C C C c C C C C c C C C C c C C C C do i=1,nn point(i)=point(i)+dstep(i) enddo kette(jnr+1,1)=point(1)*btoa kette(jnr+1,2)=point(2)*btoa kette(jnr+1,3)=point(3)/torad rewind 32 write(32,*) (kette(jnr+1,i),i=1,nn) cccccccccccccccccccccccccccccccccccccccccccccccccc jnr=jnr+1 GOTO 55 else c goto vrisearch in part 3, or stop in part 5 istatus=0 rewind 18 write(18,*) istatus icalc=0 rewind 19 write(19,*) icalc,sgn if(iter.eq.0) then iter=1 rewind 17 write(17,*) iter,iterNT DO 891 k=1,le+1 do i=1,nn rkette(k,i)=kette(le+2-k,i) enddo 891 continue DO 892 k=1,le+1 do i=1,nn kette(k,i)=rkette(k,i) enddo 892 continue rewind 25 write(25,2)(-1.0d0*tang(i),i=1,nn) rewind 26 write(26,2)(-1.0d0*t1(i),i=1,nn) endif goto 66 ENDIF ccc Abbruch -- condition end c ENDIF ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccCCCCCCCCCCCCCCCCCC tt=dmnorm(dstep,gmat,nn) write(44,*) ' dstep is ' ccccc Notbremse -- emergency break : factor=1.0d0 if(tt .gt. 0.05d0) factor=0.05d0/tt write(44,2) (dstep(I)*factor*stl,I=1,nn) DO 20 i=1,nn point(i)=point(i)+dstep(i)*factor*stl 20 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC WRITE(44,*) 'point nach dstep addition' WRITE(44,*) (point(I),I=1,nn) kette(jnr,1)=point(1)*btoa kette(jnr,2)=point(2)*btoa kette(jnr,3)=point(3)/torad write(44,*) 'Nr j in vri2 is =',jnr write(*,*) 'j vri2 = ',jnr,' iterNT = ',iterNT write(*,*) ' |dstep| ', tt,' rg ',rgnorm write(44,*) ' |dstep| ', tt,' rg ',rgnorm iterNT=iterNT+1 rewind 17 write(17,*) iter,iterNT ccc 55 continue ccc rewind 9 WRITE(9,*) jnr ccc 66 continue ccc 202 FORMAT(I5,3X,3F20.14) if(icalc.eq.1) then DO 88 k=1,le+1 88 WRITE(44,202) k, (kette(k,i),i=1,nn) endif rewind 8 DO 89 k=1,le+1 89 WRITE(8,2) (kette(k,i),i=1,nn) rewind 32 WRITE(32,2) (kette(jnr,i),i=1,nn) if(jnr.eq.le+1 .and. iterNT.eq.1 ) then c make Mma input file !! CH and CN coords interchanged WRITE(7,*)'{ (* iter = ',iter,' *) ' DO 606 k=1,le 606 WRITE(7,632) kette(k,2),kette(k,1),kette(k,3) k=le+1 WRITE(7,633) kette(k,2),kette(k,1),kette(k,3) WRITE(7,*)'}' 632 FORMAT('{',F20.14,' , ',F20.14,' , ',F20.14,'/90.0 },') 633 FORMAT('{',F20.14,' , ',F20.14,' , ',F20.14,'/90.0 }' ) endif rewind 18 write(18,*) istatus rewind 19 write(19,*) icalc,sgn call exit(istatus) END ! ! Legt gerade "kette" der Laenge "le" zwischen die Punkte "x" und "y". ! SUBROUTINE straight_chain(kette,x,y,nn,le) DOUBLE PRECISION kette(le+1,nn),x(nn),y(nn) DO 20 j=1,le+1 DO 10 i=1,nn kette(j,i)=x(i)+((j-1)*(y(i)-x(i)))/le 10 CONTINUE 20 CONTINUE END cccccccccccccccccccccccccccccccccccccccccccccccc ! Berechnet den Projektionsoperator "pr" ! aus der Suchrichtung "r" covariant, und "rc" contravariat mit ! pr= (Id-r*rc^T), pr is a mixed tensor SUBROUTINE projector(pr,r,rc,nn) DOUBLE PRECISION r(nn),rc(nn),pr(nn,nn) DOUBLE PRECISION nrm,vec_ddot c nrm=1.d0*norm(r,nn) nrm=dsqrt(vec_ddot(r,rc,nn)) DO 5 i=1,nn r(i)=r(i)/nrm rc(i)=rc(i)/nrm 5 CONTINUE DO 20 i=1,nn DO 10 j=1,nn IF (i.EQ.j) THEN pr(i,j)=-r(i)*rc(j)+1.d0 ELSE pr(i,j)=-r(i)*rc(j) ENDIF 10 CONTINUE 20 CONTINUE END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine ginvers(in,NAT,out,omat_inv) c calculates the inverse matrix of IN real*8 in,out logical omat_inv integer NAT,i,j,k dimension in(NAT,NAT),out(NAT,NAT) do 1,k=1,NAT do 1,j=1,NAT 1 out(j,k)=in(j,k) c do 10,k=1,NAT if (out(k,k).eq.0.0d0) then omat_inv=.false. return endif out(k,k)=1.0d0/out(k,k) do 20,i=1,NAT if (i.ne.k) out(i,k)=-(out(i,k)*out(k,k)) 20 continue do 30,j=1,NAT if (j.ne.k) then do 40,i=1,NAT if (i.ne.k) out(i,j)=out(i,j)+out(i,k)*out(k,j) 40 continue out(k,j)=out(k,j)*out(k,k) endif 30 continue 10 continue omat_inv=.true. return end Ccccccccccccccccccccccccccccccccccccccccccccc ! Berechnet adjunkte Matrix A zu M for n=3 SUBROUTINE adjoint(A,M,nn) DOUBLE PRECISION A(nn,nn),M(nn,nn) A(1,1)=M(2,2)*M(3,3)-M(2,3)*M(3,2) A(1,2)=M(1,3)*M(3,2)-M(1,2)*M(3,3) A(1,3)=M(1,2)*M(2,3)-M(1,3)*M(2,2) A(2,1)=M(2,3)*M(3,1)-M(2,1)*M(3,3) A(2,2)=M(1,1)*M(3,3)-M(1,3)*M(3,1) A(2,3)=M(1,3)*M(2,1)-M(1,1)*M(2,3) A(3,1)=M(2,1)*M(3,2)-M(2,2)*M(3,1) A(3,2)=M(1,2)*M(3,1)-M(1,1)*M(3,2) A(3,3)=M(1,1)*M(2,2)-M(1,2)*M(2,1) c case n=2 c A(1,1)=M(2,2) c A(1,2)=-M(1,2) c A(2,1)=-M(2,1) c A(2,2)=M(1,1) 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 cccccccccccccccccccccccccccccccc ! Multiplikation Ax=b SUBROUTINE multiply(b,A,x,nn) DOUBLE PRECISION b(nn),A(nn,nn),x(nn) DO 20 i=1,nn b(i)=0 DO 10 j=1,nn b(i)=b(i)+A(i,j)*x(j) 10 CONTINUE 20 CONTINUE END ccccccccccccccccccccccccccccccccccccccc ! Gibt die Norm von v aus DOUBLE PRECISION FUNCTION norm(v,nn) DOUBLE PRECISION v(nn),s s=0.d0 DO 10 j=1,nn s=s+v(j)**2 10 CONTINUE IF (s.LT.1.D-23) THEN norm=0.d0 ELSE norm=DSQRT(s) ENDIF END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 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 c note: dmgs is the Gram-Schmidt method 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 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 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 ccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine gettang(curveq,gmat,tangent,nn) c This routine calculate the tangent on a curve c whereas the tangent is unique determined by: c H*t=0 c |t|=1 c det[H,t]>0 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) 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) call vec_dmult_constant + (tangent,nn,1.0d0/dmnorm(tangent,gmat,nn),tangent) do 10,j = 1,nn-1 do 20,i = 1,nn ss1(j,i) = curveq(j,i) 20 continue 10 continue do 30,i = 1,nn 30 ss1(nn,i)=tangent(i) if (det(ss1,nn).lt.zero) then call vec_dmult_constant(tangent,nn,-1.0d0,tangent) endif return end ccccccccccccccccccccccccccccccccccccccccccccccccccccc 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 cccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine qr(gex,q,n) C QR-Zerlegung von 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 write(44,*) 'in QR ' do i=1,n+1 write(44,*) (q(i,j),j=1,n+1) enddo return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccc 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 ccccccccccccccccccccccccccccccccccccccccccccccccccccc 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 ccccccccccccccccccccccccccccccccccccccccccccccccccc 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 ccccccccccccccccccccccccccccccccccccccccccc subroutine buildk(dk,req,tang,rgr,nn) c build the K-Matrix real*8 dk,req,tang,rgr 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 dk(nn,nn+1)=0.0d0 return end cccccccccccccccccccccccccccccccccccccc 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) write(44,*) 'in linsolve ' do k=1,n write(44,*) (dre(k,i),i=1,n+1) enddo call loesung(dre,x,n) return end cccccccccccccccccccccccccccccccccccccccccccc 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 cccccccccccccccccccccccccccccccccccccccccccccccccc 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 ccccccccccccccccccccccccccccccccccccc 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 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 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 cccccccccccccccccccccccccccccccccccccccccccccccccccc 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 ccccccccccccccccccccccccccccccccccccc 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 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 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 ccccccccccccccccccccccccccccccccccccccccccccccc 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 cccccccccccccccccccccccccccccccccccccccccccc 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 cccccccccccccccccccccccccccccccccccccccccccccccccc 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 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 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 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 cccccccccccccccccccccccccccccccccccccccccccccccccc 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 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 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine eigen (a,m,n,d,vec,e,iff) c eigenvalues and -vectors of A real*8 a,d,vec,e,eps,tol,h,g,s,f,b,p,r,c integer m,iff,n,i,j,ni,ii,l,k,j1 dimension a(m,m), d(m), vec(m,m) dimension e(m) eps = 0.278d-16 tol = 0.211758d-21 if (n.eq.1) then d(1) = a(1,1) vec(1,1) = 1.0d0 else do 30 i = 1 , n c householder's reduction c simulation of loop do 150 i=n,2,(-1) do 20 j = 1 , i vec(i,j) = a(i,j) 20 continue 30 continue do 120 ni = 2 , n ii = n + 2 - ni do 110 i = ii , ii l = i - 2 h = 0.0d0 g = vec(i,i-1) if (l.gt.0) then do 40 k = 1 , l h = h + vec(i,k)**2 40 continue s = h + g*g if (s.lt.tol) then h = 0.0d0 else if (h.gt.0) then l = l + 1 f = g g = dsqrt(s) if (f.gt.0) then g = -g end if h = s - f*g vec(i,i-1) = f - g f = 0.0d0 do 70 j = 1 , l vec(j,i) = vec(i,j)/h s = 0.0d0 do 50 k = 1 , j s = s + vec(j,k)*vec(i,k) 50 continue j1 = j + 1 if (j1.le.l) then do 60 k = j1 , l s = s + vec(k,j)*vec(i,k) 60 continue end if e(j) = s/h f = f + s*vec(j,i) 70 continue f = f/(h+h) do 80 j = 1 , l e(j) = e(j) - f*vec(i,j) 80 continue do 100 j = 1 , l f = vec(i,j) s = e(j) do 90 k = 1 , j vec(j,k) = vec(j,k) - f*e(k) - vec(i,k)*s 90 continue 100 continue end if end if c accumulation of transformation matrices d(i) = h e(i-1) = g 110 continue 120 continue d(1) = vec(1,1) vec(1,1) = 1.0d0 do 170 i = 2 , n l = i - 1 if (d(i).gt.0) then do 150 j = 1 , l s = 0.0d0 do 130 k = 1 , l s = s + vec(i,k)*vec(k,j) 130 continue do 140 k = 1 , l vec(k,j) = vec(k,j) - s*vec(k,i) 140 continue 150 continue end if d(i) = vec(i,i) vec(i,i) = 1.0d0 do 160 j = 1 , l c diagonalization of the tridiagonal matrix vec(i,j) = 0.0d0 vec(j,i) = 0.0d0 160 continue 170 continue b = 0.0d0 f = 0.0d0 e(n) = 0.0d0 do 250 l = 1 , n c test for splitting h = eps*(dabs(d(l))+dabs(e(l))) if (h.gt.b) b = h do 180 j = l , n c test for convergence if (dabs(e(j)).le.b) go to 190 180 continue 190 if (j.eq.l) then d(l) = d(l) + f else 200 p = (d(l+1)-d(l))*0.5d0/e(l) r = dsqrt(p*p+1.0d0) if (p.lt.0) then p = p - r else p = p + r end if h = d(l) - e(l)/p do 210 i = l , n c qr transformation d(i) = d(i) - h 210 continue f = f + h p = d(j) c simulation of loop do 330 i=j-1,l,(-1) c = 1.0d0 s = 0.0d0 j1 = j - 1 do 240 ni = l , j1 ii = l + j1 - ni do 230 i = ii , ii c protection against underflow of exponents g = c*e(i) h = c*p if (dabs(p).lt.dabs(e(i))) then c = p/e(i) r = dsqrt(c*c+1.0d0) e(i+1) = s*e(i)*r s = 1.0d0/r c = c/r else c = e(i)/p r = dsqrt(c*c+1.0d0) e(i+1) = s*p*r s = c/r c = 1.0d0/r end if p = c*d(i) - s*g d(i+1) = h + s*(c*g+s*d(i)) do 220 k = 1 , n h = vec(k,i+1) vec(k,i+1) = vec(k,i)*s + h*c vec(k,i) = vec(k,i)*c - h*s 220 continue 230 continue 240 continue e(l) = s*p c convergence d(l) = c*p if (dabs(e(l)).gt.b) go to 200 c ordering of eigenvalues d(l) = d(l) + f end if 250 continue if (iff.lt.1) then ni = n - 1 do 280 i = 1 , ni k = i p = d(i) j1 = i + 1 do 260 j = j1 , n if (d(j).lt.p) then k = j p = d(j) end if 260 continue if (k.ne.i) then d(k) = d(i) d(i) = p do 270 j = 1 , n p = vec(j,i) vec(j,i) = vec(j,k) vec(j,k) = p 270 continue end if c special treatment of case n = 1 280 continue end if end if return end ccccccccccccccccccccccccccccccccccccccccccc subroutine mat_adj(matin,ndim,matout) c adjoint matrix of matin real*8 matin,matout,matscr,det integer ndim,i,j,k,l,ii,jj dimension matin(ndim,ndim),matout(ndim,ndim), + matscr(ndim-1,ndim-1) call vec_dinit(matscr,(ndim-1)*(ndim-1),0.0d0) do 10, i=1 , ndim do 20, j=1 , ndim do 30, ii=1 , ndim-1 do 40, jj=1 , ndim-1 if (ii.lt.i) then k=ii else k=ii+1 endif if (jj.lt.j) then l=jj else l=jj+1 endif matscr(ii,jj)=matin(k,l) c if(dabs(matscr(ii,jj)).lt.1.0d-30) then c matscr(ii,jj)=0.0d0 c endif 40 continue 30 continue matout(i,j)=det(matscr,ndim-1) if(dabs(matout(i,j)).lt.1.0d-30) goto 20 matout(i,j)= (-1.0d0)**(i+j)*matout(i,j) 20 continue 10 continue return end cccccccccccccccccccccccccccccccccccccccccccccccccc