PROGRAM FANTACROSS c A program to evaluate cross-correlation relaxation c INPUT PARAMETERS TO BE GIVEN AT !@!: c 3 Euler angles c 2 anisotropies of the chi tensor c 3 coordinates: position of the metal ion INCLUDE 'fantacross.h' dimension npt(MAXFE,6), * ccptx(MAXSTR,MAXFE,6), ccpty(MAXSTR,MAXFE,6), * ccptz(MAXSTR,MAXFE,6), * vx(3), vy(3),vz(3) character label*3,numero*5 common /label/ label(300,300),numero(300,300) character namepara*4 character yes*1 CHARACTER varerr*80 ! Error management string CHARACTER errmsg*80 ! Error management string logical InECHO parameter (InECHO=.true.) ! INPUT ECHO FLAG (echo/silent) logical foundAT(MAXFE,3) ! Found-reference-systems-atoms flags do i=1,MAXFE do j=1,3 foundAT(i,j)=.false. enddo enddo gltoll=0.0 print *, ' ' print *, 'This is FantaOrient' c c To be sure we have only one bloody tensor c nfe=1 c c c ******************** c *** MANUAL INPUT *** c ******************** print *, ' ' print *, 'Name of observed deltaJ file' read (*,'(a)') filename1 if(InECHO) print '(3x,a)', filename1 print *, 'Name of coordinate file' read (*,'(a)') filename2 if(InECHO) print '(3x,a)', filename2 print*, 'Input file format (1:dyana cor, 2:amber pdb)' read (*,'(i1)') nff if(InECHO) then if(nff.eq.1) then print *, 'dyana' else if(nff.eq.2) then print *, 'amber' else errmsg=' Invalid input parameter' varerr=' code:neg' goto 22 endif endif print *, 'Name for processed file' read (*,'(a)') filename3 if(InECHO) print '(3x,a)', filename3 print *, 'Name for output file' read (*,'(a)') fileout if(InECHO) print '(3x,a)', fileout print *,'Number of atoms in each structure (not more than' *,MAXDEF,'):' read (*,'(i4)') nat if(InECHO) print '(3x,i4)', nat if(nat.gt.MAXDEF) then errmsg=' Too many atoms in each structure' varerr=' code:MAXDEF' goto 22 else if (nat.lt.1) then errmsg=' Invalid input parameter' varerr=' code:neg' goto 22 endif print *, 'Number of structures (not more than',MAXSTR,'):' read (*,'(i4)') nstr if(InECHO) print '(3x,i4)', nstr if(nstr.gt.MAXSTR) then errmsg=' Too many structures required' varerr=' code:MAXSTR' goto 22 else if (nstr.lt.1) then errmsg=' Invalid input parameter' varerr=' code:neg' goto 22 endif print *, 'Value of tolerance (0-1)' read (*, '(f8.3)') PERC if(InECHO) print '(3x,f8.3)', PERC if(PERC.lt.0.or.PERC.gt.1) then errmsg=' Value of tolerance out of range' varerr=' code:range' goto 22 endif print *,' ' print *,' Would you define the reference system? (y/n) ' print *,' ' read (*,'(a1)') yes if (yes.eq.'y'.or.yes.eq.'Y') then intsys=1 nsystem=1 else intsys=0 nsystem=1 endif c select number of residue c c numero dei tre atomi npt(j,i) c if (intsys.eq.1) then do j=1,nsystem print *,'Defining reference system nr.',j,'/',nsystem,':' print *,' Insert three atom numbers:' print *,' - first fixes the x direction' read(*,'(i5)') npt(j,1) if(InECHO) print '(5x,i5)', npt(j,1) print *,' - second is the system origin' read(*,'(i5)') npt(j,2) if(InECHO) print '(5x,i5)', npt(j,2) print *,' - third is an other point of the plane x-y' read(*,'(i5)') npt(j,3) if(InECHO) print '(5x,i5)', npt(j,3) print *,' ' enddo endif c griglia sulle condizioni iniziali print *,'Number of computed solutions:' read (*,'(i4)') ngrid if(InECHO) print '(3x,i4)', ngrid c ***** READ OBSERVED DeltaJs ****** c c open (1,file=filename1,status='old') i=0 numS=0 numD=0 numT=0 print *,' ***** READ OBSERVED DeltaJs' do while (1.gt.0) read(1,'(A)',end=99) line i=i+1 read(line, * '(i3,1x,a4,1x,a4,1x,1x,f7.3,1x,f5.2,1x,f5.2)') * numres(i),namres(i),namat(i),obs(i), * tolprot(i),wprot(i) if(namat(i)(1:1).eq.' ') then namat(i)(1:3)=namat(i)(2:4) namat(i)(4:4)=' ' end if if (tolprot(i).eq.0) tolprot(i)=gltoll if ((PERC.gt.0).or.(PERC.le.1)) then ttperc=abs(obs(i)*PERC) if (ttperc.gt.tolprot(i)) then tolprot(i)=ttperc end if end if if (wprot(i).eq.0.) wprot(i)=1.0 if (wprot(i).eq.-1.) wprot(i)=0.0 enddo 99 ihp=i print *, ' ***** TOTAL No. NHs ',ihp open (2,file=filename2,status='old') open (3,file=filename3) c ***** READ COORDINATE FILE ***** c c numberatom c nameatom c nores c xa,ya,za c print *,' ***** READ COORDINATE FILE' if (nff.eq.1) then formato='(i5,1x,a4,1x,i6,6x,3f11.4)' formato1='(i5,1x,a4,1x,i6,1x,a3,2x,3f11.4)' else formato='(5x,i6,1x,a4,6x,i4,4x,3f8.3)' formato1='(5x,i6,1x,a4,1x,a3,1x,i5,4x,3f8.3)' endif ncontpoint=0 do 10 k=1,nstr do 20 i=1,nat read(2,'(a)',end=98) line read(line,formato) * numberatom,nameatom,nores,xa,ya,za c write(6,formato) c * numberatom,nameatom,nores,xa,ya,za c c c Selezione atomi per il sistema di riferimento c ccptx(k,j,l) ccpty(k,j,l) ccptz(k,j,l) c c if (intsys.eq.1) then do 15 j=1,nsystem do 15 l=1,3 if (numberatom.eq.npt(j,l)+(k-1)*nat) then ccptx(k,j,l)=xa ccpty(k,j,l)=ya ccptz(k,j,l)=za foundAT(j,l)=.true. endif 15 continue endif c c selezione degli N, HN c xn(k,i),yn(k,i),zn(k,i) c xp(k,i),yp(k,i),zp(k,i) c if (nameatom(1:1).eq.'N'.or.nameatom(2:2).eq.'N') then if(nff.eq.1) then read(line,formato1) * num_at(k,i),nam_at(k,i),num_res(k,i),nam_res(k,i), * xn(k,i),yn(k,i),zn(k,i) else read(line,formato1) * num_at(k,i),nam_at(k,i),nam_res(k,i),num_res(k,i), * xn(k,i),yn(k,i),zn(k,i) end if if(nam_at(k,i)(1:1).eq.' ') then nam_at(k,i)(1:3)=nam_at(k,i)(2:4) nam_at(k,i)(4:4)=' ' end if end if if (nameatom(1:2).eq.'HN'.or.nameatom(2:3).eq.'HN') then if(nff.eq.1) then read(line,formato1) * num_at(k,i),nam_at(k,i),num_res(k,i),nam_res(k,i), * xp(k,i),yp(k,i),zp(k,i) else read(line,formato1) * num_at(k,i),nam_at(k,i),nam_res(k,i),num_res(k,i), * xp(k,i),yp(k,i),zp(k,i) end if if(nam_at(k,i)(1:1).eq.' ') then nam_at(k,i)(1:3)=nam_at(k,i)(2:4) nam_at(k,i)(4:4)=' ' end if end if 20 continue c DATA CONTROLS: were the atoms really present? if (intsys.eq.1) then ! Reference systems were required do 16 j=1,nsystem ! Atoms defining the reference systems do 16 l=1,3 if (.not.foundAT(j,l)) then print *,' *** ERROR: The required atom',npt(j,l), * 'is not present in structure nr.',k stop endif 16 continue endif 10 continue goto 25 c Error management 22 print *,' *** ERROR: ',errmsg print *,varerr print *,' *** PROGRAM STOP.' print *,' ' STOP 25 continue c c ***** CALCULATE THE REFERENCE SYSTEM **** c c c c If an internal reference system has been selected.... c c if (intsys.eq.1) then print *,' ***** CALCULATE THE REFERENCE SYSTEM' c calcola i tre punti medi do 30 j=1,nsystem do 30 m=1,3 rccptx(j,m) = 0.0 rccpty(j,m) = 0.0 rccptz(j,m) = 0.0 30 continue do 40 j=1,nsystem do 40 k=1,nstr do 40 m=1,3 rccptx(j,m) = rccptx(j,m) + ccptx(k,j,m)/dble(nstr) rccpty(j,m) = rccpty(j,m) + ccpty(k,j,m)/dble(nstr) rccptz(j,m) = rccptz(j,m) + ccptz(k,j,m)/dble(nstr) 40 continue c c calcolo versore Z, X ed Y c do 50 j=1,nsystem vz(1)=(rccpty(j,2)-rccpty(j,1))* * (rccptz(j,3)-rccptz(j,1))- * (rccpty(j,3)-rccpty(j,1))*(rccptz(j,2)-rccptz(j,1)) vz(2)=(rccptx(j,3)-rccptx(j,1))* * (rccptz(j,2)-rccptz(j,1))- * (rccptx(j,2)-rccptx(j,1))*(rccptz(j,3)-rccptz(j,1)) vz(3)=(rccptx(j,2)-rccptx(j,1))* * (rccpty(j,3)-rccpty(j,1))- * (rccptx(j,3)-rccptx(j,1))*(rccpty(j,2)-rccpty(j,1)) rvz=sqrt(vz(1)**2+vz(2)**2+vz(3)**2) vz(1)=vz(1)/rvz vz(2)=vz(2)/rvz vz(3)=vz(3)/rvz c c calcolo versore X = PT(1) - PT(2) c e' orientato dal primo al secondo che e' il centro c vx(1)=rccptx(j,1)-rccptx(j,2) vx(2)=rccpty(j,1)-rccpty(j,2) vx(3)=rccptz(j,1)-rccptz(j,2) rvx=sqrt(vx(1)**2+vx(2)**2+vx(3)**2) vx(1)=vx(1)/rvx vx(2)=vx(2)/rvx vx(3)=vx(3)/rvx c c calcolo versore Y DESTROSO c vy(1)=vz(2)*vx(3)-vx(2)*vz(3) vy(2)=vx(1)*vz(3)-vz(1)*vx(3) vy(3)=vz(1)*vx(2)-vx(1)*vz(2) c c vettori da visualizzare c center(j,1)=rccptx(j,2) center(j,2)=rccpty(j,2) center(j,3)=rccptz(j,2) ssx(j,1)=rccptx(j,2)+vx(1) ssx(j,2)=rccpty(j,2)+vx(2) ssx(j,3)=rccptz(j,2)+vx(3) ssy(j,1)=rccptx(j,2)+vy(1) ssy(j,2)=rccpty(j,2)+vy(2) ssy(j,3)=rccptz(j,2)+vy(3) ssz(j,1)=rccptx(j,2)+vz(1) ssz(j,2)=rccpty(j,2)+vz(2) ssz(j,3)=rccptz(j,2)+vz(3) c c matrici con i coseni direttori dei sistemi di rifer. c axsys(j,1)=vx(1) axsys(j,2)=vx(2) axsys(j,3)=vx(3) aysys(j,1)=vy(1) aysys(j,2)=vy(2) aysys(j,3)=vy(3) azsys(j,1)=vz(1) azsys(j,2)=vz(2) azsys(j,3)=vz(3) 50 continue c c c Else, we put the reference in (0,0,0), i.e. in the origin of c the lab system. The numbers relative to the c orientation of the tensor will be thus c dependent on the orientation of the protein with respect to the c lab frame in the input file. c c else rccptx(1,2) = 0. rccpty(1,2) = 0. rccptz(1,2) = 0. c c matrici con i coseni direttori di sistema di rifer. c do 70 j=1,nfe axsys(j,1)=1. axsys(j,2)=0. axsys(j,3)=0. aysys(j,1)=0. aysys(j,2)=1. aysys(j,3)=0. azsys(j,1)=0. azsys(j,2)=0. azsys(j,3)=1. c c vettori da visualizzare c center(j,1)=rccptx(j,2) center(j,2)=rccpty(j,2) center(j,3)=rccptz(j,2) ssx(j,1)=rccptx(j,2)+1. ssx(j,2)=rccpty(j,2) ssx(j,3)=rccptz(j,2) ssy(j,1)=rccptx(j,2) ssy(j,2)=rccpty(j,2)+1. ssy(j,3)=rccptz(j,2) ssz(j,1)=rccptx(j,2) ssz(j,2)=rccpty(j,2) ssz(j,3)=rccptz(j,2)+1. 70 continue endif c c We're done with the reference system, so now we move on to c processing the input file. c c c ***** WRITE PROCESSED COORDINATE + DeltaJ FILE **** c print *,' ***** WRITE PROCESSED INPUT FILE' do 140 k=1,nstr icontprot=0 icont15n=0 do 150 j=1,ihp do 160 i=1,nat if (num_res(k,i).eq.numres(j)) then if (nam_at(k,i).eq.namat(j) * .or.nam_at(k,i)(2:4).eq.namat(j)(1:3) * .or.nam_at(k,i)(1:3).eq.namat(j)(2:4)) then icontprot=icontprot+1 write(3, * '(i6,1x,a4,1x,i3,1x,a3,1x,3F8.3,1x,F7.3,2(1x,f7.2))') * num_at(k,i),nam_at(k,i),num_res(k,i),nam_res(k,i), * xp(k,i),yp(k,i),zp(k,i),obs(j), * tolprot(j),wprot(j) do ii=1,nat if (num_res(k,ii).eq.numres(j).and.nam_at(k,i) * (2:3).eq.nam_at(k,ii)) then icont15n=icont15n+1 write(3, * '(i6,1x,a4,1x,i3,1x,a3,1x,3F8.3,1x,F7.3,2(1x,f7.2))') * num_at(k,ii),nam_at(k,ii),num_res(k,ii),nam_res(k,ii), * xn(k,ii),yn(k,ii),zn(k,ii),obs(j), * tolprot(j),wprot(j) goto 150 end if end do print*, '1',nam_at(k,i)(2:3),numres(j),namres(j) goto 150 end if end if 160 continue print*, '2',namat(j),numres(j),namres(j) 150 continue 140 continue 98 continue if (icontprot.lt.ihp) then print *,' ****** ERROR: SOME ATOMS NOT FOUND ' print *,' ' print *,' ****** THE PROGRAM WILL STOP' print *,' ' STOP endif if (icontprot.ne.icont15n) then print *,' ****** ERROR: No. OF 15N ATOMS NOT EQUAL', * ' TO No. OF HN' print *,' ' print *,' ****** THE PROGRAM WILL STOP' print *,' ' print *,icontprot,icont15n STOP endif close(1) close(2) close(3) nhp=icontprot*2 open (1,file=filename3,status='old') c c ***** READ PROCESSED INPUT FILE ***** c coord. dei protoni cx(k,i),cy(k,i),cz(k,i) c valori dei param. di shift obs tolprot wprot c i=0 c print *,' ***** READ PROCESSED INPUT FILE' c c nhp is twice the number of constraints c With i odd, the variables read with the instruction below c refer to an HN, whereas the values with index i+1 (thus, even) c refer to the relative 15N. Note that the values of obs, c tolprot and wprot are the same for i and i+1. c do 170 k=1,nstr do 170 i=1,nhp read(1,'(A)',end=199) line read(line,'(7x,A3,A5,5x,3F8.3,1x,F7.3,2(1x,f7.2))') * label(k,i),numero(k,i),cx(k,i),cy(k,i),cz(k,i), * obs((k-1)*nhp+i), * tolprot((k-1)*nhp+i),wprot((k-1)*nhp+i) 170 continue 199 continue close(1) open (3,file=fileout) write(3,'(a)') ' ' write(3,'(1x,a,a)') 'Input file = ',filename3 write(3,'(1x,a,i6)') 'Total protons ',nhp/2*nstr write(3,'(a)') ' ' write(3,'(a)') ' ***** SOLUTIONS *****' c c ***** LUNCH SIMPLEX ****** c print *,' ' print *,' ' print *,' RUNNING SIMPLEX..........' print *,' ' c print *,' ***** RESULTS' print *,' ' call simplexrun() end C--------------------------------------------------------------- C C SIMPLEXRUN C C Calculates the susceptibilty tensor of paramagnetic system C which minimize the experimental-to-calculated squared shift C difference. C--------------------------------------------------------------- SUBROUTINE SIMPLEXRUN() C include 'fantacross.h' DIMENSION SIMP(MP,NP),EXTR(NP),VAL(MP) c dimension axx(MAXFE,3), ayy(MAXFE,3), azz(MAXFE,3) c inizializzazione SIMPLESSO OLDRESID=1.e+9 RESID=2.e+9 iviolation=2000000 ioldvio=1000000 c griglia sulle condizioni iniziali delta=3.14159/dble(ngrid+1) do 1000 l=1,ngrid do 999 nk=1,nfe SIMP(1,(nk-1)*MPAR+1)=l*delta*(1.-2*rand()) SIMP(1,(nk-1)*MPAR+2)=l*delta*(1.-2*rand()) SIMP(1,(nk-1)*MPAR+3)=l*delta*(1.-2*rand()) simp(1,(nk-1)*MPAR+4)=1000.*(1.-2*rand()) simp(1,(nk-1)*MPAR+5)=1000.*(1.-2*rand()) c SIMP(1,1)=0 c SIMP(1,2)=0 c SIMP(1,3)=0 c SIMP(1,4)=2 c SIMP(1,5)=0.5 999 continue DO 10 I=1,MP DO 10 J=1,NP SIMP(I,J)=SIMP(1,J) if (i.eq.1) then SIMP(1,J)=SIMP(1,J)*0.8 endif if (i.eq.j+1) then SIMP(I,J)=SIMP(I,J)/0.8*1.2 endif 10 CONTINUE DO 11 I=1,MP DO 12 J=1,NP 12 EXTR(J)=SIMP(I,J) VAL(I)=CASH(EXTR) 11 CONTINUE CALL SIMPCALC(SIMP,VAL) c criterio di minimizzazione della soluzione ottima if (oldresid.ge.resid) then OLDRESID=RESID ioldvio=iviolation do 13 i=1,nfe optphi(i)=phi(i) optteta(i)=teta(i) optomega(i)=omega(i) opta1(i)=a1dip(i) opta2(i)=a2dip(i) 13 continue endif 1000 continue do 998 nk=1,nfe EXTR((nk-1)*MPAR+1)=optphi(nk) EXTR((nk-1)*MPAR+2)=optteta(nk) EXTR((nk-1)*MPAR+3)=optomega(nk) EXTR((nk-1)*MPAR+4)=opta1(nk) EXTR((nk-1)*MPAR+5)=opta2(nk) 998 continue print *,' ' WRITE(*,FMT='(1X,A)') ' ***** BEST SOLUTION *****' print *,' ' WRITE (*,20) ioldvio, OLDRESID WRITE (*,21) (optphi(i),optteta(i),optomega(i),i=1,nfe) WRITE (*,22) (opta1(i),opta2(i),i=1,nfe) WRITE(3,FMT='(1X,A)') ' ' WRITE(3,FMT='(1X,A)') ' ***** BEST SOLUTION *****' WRITE(3,FMT='(1X,A)') ' ' WRITE (3,21) (optphi(i),optteta(i),optomega(i),i=1,nfe) WRITE (3,22) (opta1(i),opta2(i),i=1,nfe) c visualizzazione dei coseni direttori c if (intsys.eq.1) then do 777 i=1,nfe P=optphi(i) T=optteta(i) O=optomega(i) asxx(i,1)=cos(P)*cos(O) asxx(i,2)=sin(P)*cos(O) asxx(i,3)=sin(O) asyy(i,1)=-cos(T)*sin(P)-sin(O)*cos(P)*sin(T) asyy(i,2)=cos(T)*cos(P)-sin(O)*sin(P)*sin(T) asyy(i,3)=sin(T)*cos(O) aszz(i,1)=sin(T)*sin(P)-sin(O)*cos(P)*cos(T) aszz(i,2)=-sin(T)*cos(P)-sin(O)*sin(P)*cos(T) aszz(i,3)=cos(T)*cos(O) 777 continue c endif CALL PRRES 20 FORMAT (1X,' VIOLATIONS ',I6,1X,' RESIDUALS ',F10.3) 21 FORMAT(1X,'PHI= ',F9.3,1X,'THETA= ',F9.3,1X,'OMEGA= ',F9.3) 22 FORMAT (1X,'A1= ',F9.3,1X,'A2= ',F9.3) ERR=CASH(EXTR) call DISPSHIFT RETURN END C------------------------------------------------------------- C C SIMPCALC C C-------------------------------------------------------------- SUBROUTINE SIMPCALC(P,Y) C include 'fantacross.h' PARAMETER (AL=1.,BE=.5,GA=2.) DIMENSION P(MP,NP),Y(MP),PR(MP),PT(MP),PE(MP), * PC(MP),PK(MP),PBAR(MP),OLDP(NP),EXTR(MP) c inizializzazione dell'algoritmo ITER=0 RNDOM=0 1 perm=1 do 1000 while (perm.gt.0) perm=0 do 1100 i=1,np if (Y(i).gt.Y(i+1)) then app=Y(i+1) Y(i+1)=Y(i) Y(i)=app do 1200 j=1,np app=P(i+1,j) P(i+1,j)=P(i,j) P(i,j)=app 1200 continue perm=perm+1 endif 1100 continue 1000 continue AVERAGE=0 ASTDEV=0 DO 100 I=1,MP 100 AVERAGE=AVERAGE+Y(I) AVERAGE=AVERAGE/MP DO 101 I=1, MP 101 ASTDEV=ASTDEV+(Y(I)-AVERAGE)**2. ERR=DSQRT(ASTDEV/MP) RESID=Y(1) IF (ITER.gt.39999.OR.ERR.lt.1.D-9) THEN WRITE (*,200) ITER,IVIOLATION,RESID WRITE (3,200) ITER,IVIOLATION,RESID WRITE (*,201) * (P(1,(i-1)*MPAR+1), * P(1,(i-1)*MPAR+2), * P(1,(i-1)*MPAR+3), * P(1,(i-1)*MPAR+4),P(1,(i-1)*MPAR+5),i=1,nfe) WRITE (3,201) * (P(1,(i-1)*MPAR+1), * P(1,(i-1)*MPAR+2), * P(1,(i-1)*MPAR+3), * P(1,(i-1)*MPAR+4),P(1,(i-1)*MPAR+5),i=1,nfe) 200 FORMAT (1X,'ITERS ',I6,1X,' VIOLATIONS ',I6, * 1X,'RESIDUALS= ',F10.3) 201 FORMAT(1X,'PHI= ',F9.3,1X,'THETA= ',F9.3,1X,'OMEGA= ',F9.3, * 1X,'A1= ',F9.3,1X,'A2= ',F9.3) do 997 nk=1,nfe phi(nk)=P(1,(nk-1)*MPAR+1) teta(nk)=P(1,(nk-1)*MPAR+2) omega(nk)=P(1,(nk-1)*MPAR+3) a1dip(nk)=P(1,(nk-1)*MPAR+4) a2dip(nk)=P(1,(nk-1)*MPAR+5) 997 continue RETURN ENDIF C-------------------------------------------------------- ITER=ITER+1 DO 12 J=1,NP 12 PBAR(J)=0 DO 13 I=1,NP DO 13 J=1,NP 13 PBAR(J)=PBAR(J)+P(I,J) DO 14 J=1,NP PBAR(J)=PBAR(J)/NP PR(J)=(1.+AL)*PBAR(J)-AL*P(MP,J) 14 CONTINUE YR=CASH(PR) DO 15 J=1,NP PK(J)=PR(J) 15 CONTINUE YK=YR IF(YR.lt.Y(NP)) THEN IF(YR.lt.Y(1)) THEN DO 16 J=1,NP PE(J)=GA*PR(J)+(1.-GA)*PBAR(J) 16 CONTINUE YE=CASH(PE) iter=iter+1 IF (YE.LT.Y(1)) THEN DO 17 J=1,NP PK(J)=PE(J) 17 CONTINUE YK=YE c write(*,FMT='(BN,A)') 'expand' ENDIF ENDIF ELSE DO 18 j=1,NP PT(j)=P(MP,j) 18 CONTINUE YT=Y(MP) IF (YR.LT.YT) THEN DO 19 J=1,NP PT(J)=PR(J) 19 CONTINUE YT=YR ENDIF DO 20 J=1,NP PC(J)=BE*PT(J)+(1.-BE)*PBAR(J) 20 CONTINUE YC=CASH(PC) iter=iter+1 IF (YC.lt.Y(NP)) THEN DO 21 J=1,NP PK(J)=PC(J) 21 CONTINUE YK=YC c write(*,FMT='(BN,A)') 'contract' ELSE do 22 i=1,NP OLDP(i)=P(2,i) 22 continue DO 24 I=2,NP DO 23 J=1,NP P(I,J)=.5*(P(1,J)+P(I,J)) 23 CONTINUE PR(I)=P(I,I) Y(I)=CASH(PR) 24 CONTINUE iter=iter+NP-1 do 25 i=1,NP OLDP(i)=abs(OLDP(i)-P(2,i)) 25 continue oldpoppa=0 do 251 i=1,NP-1 oldpoppa=oldpoppa+oldp(i) 251 continue if (abs(oldpoppa).lt.1.D-8.and.RNDOM.lt.10) then c iter=0 RNDOM=RNDOM+1 write(*,FMT='(BN,A,F16.6)') 'random ',Y(1) do 2001 il=1,MP do 2002 im=1,NP rmr=1+RAND() P(il,im)=P(il,im)*rmr 2002 continue 2001 continue DO 2011 Il=1,MP DO 2012 Jm=1,NP 2012 EXTR(Jm)=P(Il,Jm) Y(Il)=CASH(EXTR) 2011 CONTINUE goto 1 endif if (abs(oldpoppa).lt.1.D-9) then iter=50000 endif DO 26 J=1,NP PK(j)=0.5*(P(1,j)+P(mp,j)) 26 CONTINUE YK=CASH(PK) iter=iter+1 c write(*,FMT='(BN,A)') 'shrink' ENDIF ENDIF DO 27 J=1,NP P(MP,J)=PK(J) 27 CONTINUE Y(MP)=YK GO TO 1 END C------------------------------------------------------------- C C CASH C C-------------------------------------------------------------- FUNCTION CASH(vett) C include 'fantacross.h' character label*3,numero*5 common /label/ label(300,300),numero(300,300) PARAMETER(MBRANC=90) DIMENSION WR( MBRANC ) DIMENSION ZRR(MBRANC,MBRANC),ZRI(MBRANC,MBRANC) DIMENSION WK1(MBRANC),WK2(MBRANC),WK3(MBRANC) dimension vett(NP) dimension tx(MBRANC,MBRANC) dimension tx1(MBRANC,MBRANC) dimension txi(MBRANC,MBRANC) c calcolo dell'errore e delle violazioni IVIOLATION=0 TMP1=0 i=1 do while (i.le.nhp*nstr) shift(i)=0.0 i=i+1 enddo gj=6./7. !7./6 !4./3. !6./7. spin=2.5 !6 !15./2. !2.5 akt= 1.3807e-23*298 chiave=4*3.1428e-7*9.2741e-24**2*(gj)**2* & spin*(spin+1.)/(3*akt) ! chiave=4.6e-32 write(6,*)chiave do 1 n=1,nstr P=vett(1) T=vett(2) O=vett(3) A1D=vett(4) A2D=vett(5) ! INPUT PARAMETERS (!@!) p=1.779 !Euler angles t=-2.718 o=.787 a1d=2.076e-32 !susceptibility anisotropies a2d=-.711e-32 fmex=-5.453 !position of the metal ion fmey=18.965 fmez=9.824 !!!!!!!!!!!!!!!!!!!!!!!! chiz=2./3.*a1d+chiave chix=(a2d-chiz+3*chiave)/2 chiy=chix-a2d ! chix=chiave ! chiy=chiave ! chiz=chiave write(6,*)'chi=',chix,chiy,chiz axx=cos(P)*cos(O) axy=sin(P)*cos(O) axz=sin(O) ayx=-cos(T)*sin(P)-sin(O)*cos(P)*sin(T) ayy=cos(T)*cos(P)-sin(O)*sin(P)*sin(T) ayz=sin(T)*cos(O) azx=sin(T)*sin(P)-sin(O)*cos(P)*cos(T) azy=-sin(T)*cos(P)-sin(O)*sin(P)*cos(T) azz=cos(T)*cos(O) tmp2=0.0 i=1 ihp=(n-1)*nhp c c loop on all constraints c c We go up to nhp-1, and increment i by two c c do 10 while (i.le.nhp-1) c c First of all, we calculate the vector N-H c c c n indicates the n-th structure, i the i-th constraint c xNH=CX(n,i)-CX(n,i+1) yNH=CY(n,i)-CY(n,i+1) zNH=CZ(n,i)-CZ(n,i+1) rNH=sqrt(xNH**2+yNH**2+zNH**2) szNH=(xNH*azx+yNH*azy+zNH*azz) sxNH=(xNH*axx+yNH*axy+zNH*axz) syNH=(xNH*ayx+yNH*ayy+zNH*ayz) !C check che torna uguale facendo prima la rotazione ! write(6,*)sxNH,syNH,szNH ! cx(n,i)=cx(n,i)-fmex ! cy(n,i)=cy(n,i)-fmey ! cz(n,i)=cz(n,i)-fmez ! cx(n,i+1)=cx(n,i+1)-fmex ! cy(n,i+1)=cy(n,i+1)-fmey ! cz(n,i+1)=cz(n,i+1)-fmez ! szNH1=(cx(n,i)*azx+cy(n,i)*azy+cz(n,i)*azz) ! sxNH1=(cx(n,i)*axx+cy(n,i)*axy+cz(n,i)*axz) ! syNH1=(cx(n,i)*ayx+cy(n,i)*ayy+cz(n,i)*ayz) ! szNH2=(cx(n,i+1)*azx+cy(n,i+1)*azy+cz(n,i+1)*azz) ! sxNH2=(cx(n,i+1)*axx+cy(n,i+1)*axy+cz(n,i+1)*axz) ! syNH2=(cx(n,i+1)*ayx+cy(n,i+1)*ayy+cz(n,i+1)*ayz) ! sxNH=sxNH1-sxNH2 ! syNH=syNH1-syNH2 ! szNH=szNH1-szNH2 ! write(6,*)sxNH,syNH,szNH ! C fine check c c r should be constant, and equal to 1.02 A. c We'll use this as an internal check c if(rNH.gt.1.08.or.rNH.lt.0.94) then ! print*, '***** ERROR: r(H-N) has a wrong value: ',r ! print*, '***** Check data relative to constraint ', ! * (i+1)/2 ! print*, '***** The program should stop' ! STOP 1 end if c c Go on with calculations c xapp=CX(n,i)-fmex yapp=CY(n,i)-fmey zapp=CZ(n,i)-fmez r=sqrt(xapp**2+yapp**2+zapp**2) scalz=(xapp*azx+yapp*azy+zapp*azz) scalx=(xapp*axx+yapp*axy+zapp*axz) scaly=(xapp*ayx+yapp*ayy+zapp*ayz) tx1(1,1)=(3*scalx**2-r**2)*chix ! tx1(1,2)=3*scalx*scaly*chiy ! tx1(1,3)=3*scalx*scalz*chiz ! tx1(2,1)=3*scalx*scaly*chix tx1(2,2)=(3*scaly**2-r**2)*chiy ! tx1(2,3)=3*scaly*scalz*chiz ! tx1(3,1)=3*scalx*scalz*chix ! tx1(3,2)=3*scaly*scalz*chiy tx1(3,3)=(3*scalz**2-r**2)*chiz tave=(tx1(1,1)+tx1(2,2)+tx1(3,3))/3. tx(1,1)=tx1(1,1)-tave tx(1,2)=3*scalx*scaly*(chiy+chix)/2. tx(1,3)=3*scalx*scalz*(chiz+chix)/2. tx(2,1)=tx(1,2) tx(2,2)=tx1(2,2)-tave tx(2,3)=3*scaly*scalz*(chiz+chiy)/2. tx(3,1)=tx(1,3) tx(3,2)=tx(2,3) tx(3,3)=tx1(3,3)-tave DO 45 IJ=1,3 DO 45 J=1,3 tx(ij,j)=tx(ij,j)/(4*3.1416*r**5) ! DEVE ESSERE R**5 45 CONTINUE CALL F02AXF(TX,MBRANC,TXI,MBRANC,3,WR,ZRR,MBRANC,ZRI,MBRANC $ ,WK1,WK2,WK3,0) ! DO 46 IJ=1,3 ! write(6,*)wr(ij) ! DO 46 J=1,3 ! write(6,*)ZRR(IJ,J) !46 CONTINUE ! write(6,*)wr(1),wr(2),wr(3) ! stop gammai=2.67e8 gamman=2.71e7 gammae=1.76e11 hbar=1.055e-34 amie=9.274e-24 b0=800./42 tauc=4e-9 r6=(1e-10)**3 *rNH**3 cost=(1e-10)**3 !(per avere shift in ppm: *1e6) ! write(6,*)wr(1)/cost,wr(2)/cost,wr(3)/cost cost1=gammai**2*gamman*gammae*hbar**2/r6*2/5.*1e-7**2 cost2=gj*amie*b0/(3.*akt)*spin*(spin+1) cost3=tauc*(4.+3./(1.+(b0*42e6*2*3.1416*tauc)**2)) costnew=1.e-7/(15*3.1416)*b0*gammai**2*gamman*hbar/r6 ! & *gammae*hbar/gj/amie costetaxz=(sxNH*zrr(1,1)+syNH*zrr(2,1)+szNH*zrr(3,1))/rNH costetayz=(sxNH*zrr(1,2)+syNH*zrr(2,2)+szNH*zrr(3,2))/rNH costetazz=(sxNH*zrr(1,3)+syNH*zrr(2,3)+szNH*zrr(3,3))/rNH !potrebbero essere scambiate le righe con le colonne di zrr(i,j) anisot=(wr(1)*(3*costetaxz**2-1)+wr(2)*(3*costetayz**2-1) & +wr(3)*(3*costetazz**2-1)) write(6,130)label(n,i),numero(n,i),' CCR:', & -anisot/cost*costnew*cost3 130 format(A3,1x,A5,a15,f11.4) !controllo shift, levando -tave in matrice ! ! write(6,*)(wr(1)+wr(2)+wr(3))/3./(1e-10)**3*1e6 ! ! a1d=(chiz-(chix+chiy)/2.)/3.77e-35 a2d=(chix-chiy)/3.77e-35 g1=(sqrt(3.)*scalz-r)*(sqrt(3.)*scalz+r) g2=(scalx-scaly)*(scalx+scaly) shift(ihp+i)=shift(ihp+i)+ * (A1D*G1+1.5*A2D*G2)/r**5 i=i+2 10 continue c c Now, let's see whether the calculated DeltaJ's c agree with the input data c write(6,*)'.......... PCS ............' do 3 i=1,nhp-1,2 tmp2=abs(shift(ihp+i)-obs(ihp+i)) * -tolprot(ihp+i) if (tmp2.gt.0.0) then IVIOLATION=IVIOLATION+1 TMP1=tmp1+tmp2**2*wprot(ihp+i) endif write(6,*)label(n,i),numero(n,i),shift(ihp+i) 3 continue 1 continue stop CASH=tmp1 RETURN END C-------------------------------------------------------------- C C DISPSHIFT C C-------------------------------------------------------------- SUBROUTINE DISPSHIFT C include 'fantacross.h' character filestr*21, ss*40 dimension VV(MAXSTR,MAXOS),VMAX(maxos),VMIN(maxos) dimension nnmax(maxos),nnmin(maxos) c visualizzazione delle violazioni open (1,file=filename3,status='old') write(3,'(a)') ' ' write(3,'(a)') * '**********************************************************' write(3,'(1x,a)') * ' PROTONS MIN ERR. MAX ERR. OBS. ' rm=0 tmp1=0 smed=0 stdev=0 do 2 i=1, nstr do 2 j=1, nhp-1,2 tmp2=abs(obs(j+(i-1)*nhp)-SHIFT(j+(i-1)*nhp)) * -tolprot(j+(i-1)*nhp) if (tmp2.gt.0.0) then if (SHIFT(j+(i-1)*nhp).gt.obs(j+(i-1)*nhp)) then VV(i,j)= tmp2 else VV(i,j)= -tmp2 endif else VV(i,j)= 0. endif 2 continue do 44 i=1, nhp-1,2 VMAX(i)=0. VMIN(i)=1000. 44 continue do 4 j=1, nhp-1,2 do 4 i=1, nstr if (abs(VV(i,j)).gt.abs(VMAX(j))) then VMAX(j)=VV(i,j) c nnmax=i endif if ((abs(VV(i,j)).lt.abs(VMIN(j))).and.(VV(i,j).ne.0.)) then VMIN(j)=VV(i,j) c nnmin=i endif 4 continue do 7 j=1, nhp-1,2 do 8 i=1, nstr if (abs(VV(i,j)).gt.0.) then ss(i:i+1)="*" else ss(i:i+1)=" " endif 8 continue if (VMIN(j).eq.1000.) then VMIN(j)=0. endif read(1,'(A)') filestr c read(1,'(\)') write(3,'(1X,A,1X,F7.3,1X,F7.3,1X,F7.3,1x,A)') filestr, * VMIN(j), VMAX(j), obs(J), ss 7 continue write(3,'(1x)') write(3,'(1x)') write(3,'(1x)') write(3,'(1x,a)') * ' PROTONS MIS. CALC. ERR.^2' rewind(1) do 3 j=1,nhp*nstr-1,2 tmp2=abs(obs(j)-SHIFT(J))-tolprot(J) read(1,'(A)') filestr c read(1,'(\)') if (tmp2.gt.0.0) then IVIOLATION=IVIOLATION+1 rm=rm+1 smed=smed+tmp2**2*wprot(j) tmp1=tmp2**2*wprot(j) write(3,'(1X,A,1X,F7.3,1X,F7.3,1X,F10.5)') filestr, * obs(J),shift(j),tmp1 15 FORMAT(1X,A,1X,F7.3,1X,F7.3,1X,F7.3) else write(3,'(1X,A,1X,F7.3,1X,F7.3)') filestr, * obs(J),shift(j) endif 3 continue do 41 j=1,nhp*nstr-1,2 tmp2=abs(obs(j)-SHIFT(J))-tolprot(J) if (tmp2.gt.0.0) then stdev=stdev+(tmp2**2-med)**2*wprot(j) endif 41 continue if(rm.gt.EPS) then smed=smed/rm stdev=stdev/rm else ! No errors (this control is to avoid a 0/0) smed=0. stdev=0. endif print *,' ' WRITE(3,'(A)') ' ' write(*,11) 'mean value of error =',smed 11 format(1X,A,1X,F10.6) write(3,12) 'mean value of error =',smed 12 format(1X,A,1X,F10.6) write(*,13) 'Standard deviation =',stdev 13 format(1X,A,1X,F10.6) write(3,14) 'Standard deviation =',stdev 14 format(1X,A,1X,F10.6) close(1) close(3) RETURN END C ---------------------------------- C-------------------------------------------------------------- C C PRES C C-------------------------------------------------------------- SUBROUTINE PRRES C include 'fantacross.h' c c saxx sayy sazz : cos. dir. delle matr. di rot. c c saxsys saysys sazsys : cos. dir. dei sist. di rifer. c dimension saxx(3),sayy(3),sazz(3), * appx(3),appy(3),appz(3), * saxsys(3),saysys(3),sazsys(3) character namatom*4,namress*3 PI=3.1415927 c c open pdb file c open (2,file=filename2,status='old') i=0 do while (i.le.nstr*nat) i=i+1 read(2,'(A)',END=99) line if (nff.eq.1) then read(line,formato1) * numberatom,namatom,nores,namress,xa,ya,za else read(line,formato1) * numberatom,namatom,namress,nores,xa,ya,za endif enddo 99 i=0 A1=opta1(1) A2=opta2(1) if (nsystem.gt.1) then do 1 i=1,3 saxsys(i)=axsys(1,i) saysys(i)=aysys(1,i) sazsys(i)=azsys(1,i) 1 continue else do 2 i=1,3 saxsys(i)=axsys(1,i) saysys(i)=aysys(1,i) sazsys(i)=azsys(1,i) 2 continue endif do 3 i=1,3 saxx(i)=asxx(1,i) sayy(i)=asyy(1,i) sazz(i)=aszz(1,i) 3 continue sdzx=saxx(2)*sayy(3)-saxx(3)*sayy(2) if (sdzx/sazz(1).gt.0.0) then c destrorso print *,' ' else c sinistrorso print *,' ' do 4 i=1,3 appx(i)=saxx(i) appy(i)=sayy(i) 4 continue do 5 i=1,3 saxx(i)=appy(i) sayy(i)=appx(i) 5 continue A2=-A2 endif print *, *'***********************************************************' WRITE(3,'(A)') ' ' WRITE(3,'(A)') * '**********************************************************' numberatom=numberatom+1 namatom=' CEN' namress='PMC' nores=nores+1 c c The following line should be avoided, and all n's below c replaced by 1, because in this system we have always one tensor c but I'm too lazy... c n=1 WRITE(*,'(A)') * ' ANGLES BETWEEN REFERENCE AXIS AND TENSOR AXIS (degrees) ' WRITE(3,'(A)') ' ' WRITE(3,'(A)') * ' ANGLES BETWEEN REFERENCE AXIS AND TENSOR AXIS (degrees) ' csa=saxx(1)*saxsys(1)+saxx(2)*saxsys(2)+ * saxx(3)*saxsys(3) csb=sayy(1)*saysys(1)+sayy(2)*saysys(2)+ * sayy(3)*saysys(3) csg=sazz(1)*sazsys(1)+sazz(2)*sazsys(2)+ * sazz(3)*sazsys(3) write (*,23) acos(csa)*180./PI,acos(csb)*180./PI,acos(csg)*180./PI write (3,23) acos(csa)*180./PI,acos(csb)*180./PI,acos(csg)*180./PI WRITE(*,'(A)') ' FIVE TENSOR PARAMETERS ' WRITE(3,'(A)') ' ' WRITE(3,'(A)') ' FIVE TENSOR PARAMETERS ' CALL PANG(0,1) fattoreconv=(18.7**2-11.7**2)/(15*1.3807e-23*300)*(2.675e8*2.71e7* & 6.626e-34)/(4*3.1415**2)/(4*3.1415)*1e30 WRITE (*,22) * A1,'(',A1/fattoreconv*1e32,' m^3 e+32) ', * A2,'(',A2/fattoreconv*1e32,' m^3 e+32) ' WRITE (3,22) * A1,'(',A1/fattoreconv*1e32,' m^3 e+32) ', * A2,'(',A2/fattoreconv*1e32,' m^3 e+32) ' c c vettori da visualizzare c tsx(n,1)=rccptx(n,2)+saxx(1) tsx(n,2)=rccpty(n,2)+saxx(2) tsx(n,3)=rccptz(n,2)+saxx(3) tsy(n,1)=rccptx(n,2)+sayy(1) tsy(n,2)=rccpty(n,2)+sayy(2) tsy(n,3)=rccptz(n,2)+sayy(3) tsz(n,1)=rccptx(n,2)+sazz(1) tsz(n,2)=rccpty(n,2)+sazz(2) tsz(n,3)=rccptz(n,2)+sazz(3) c c append atoms to pdb file c numberatom=numberatom+1 namatom=' TEX' WRITE(3,'(A)') ' ' WRITE(3,'(A)') ' PDB COORDINATES OF TENSOR SYSTEM' WRITE(3,'(A4,1x,i6,1x,a4,1x,a3,2x,i4,4x,3f8.3)') * 'ATOM',numberatom,namatom,namress,nores, * tsx(1,1),tsx(1,2),tsx(1,3) numberatom=numberatom+1 namatom=' TEY' WRITE(3,'(A4,1x,i6,1x,a4,1x,a3,2x,i4,4x,3f8.3)') * 'ATOM',numberatom,namatom,namress,nores, * tsy(1,1),tsy(1,2),tsy(1,3) numberatom=numberatom+1 namatom=' TEZ' WRITE(3,'(A4,1x,i6,1x,a4,1x,a3,2x,i4,4x,3f8.3)') * 'ATOM',numberatom,namatom,namress,nores, * tsz(1,1),tsz(1,2),tsz(1,3) numberatom=numberatom-3 WRITE(*,'(A,A)') * ' ANGLES BETWEEN REFERENCE AXIS AND TENSOR AXIS ' , * '(degrees) 1 PERMUTATION' WRITE(3,'(A)') ' ' WRITE(3,'(A,A)') * ' ANGLES BETWEEN REFERENCE AXIS AND TENSOR AXIS ' , * '(degrees) 1 PERMUTATION' do 6 i=1,3 appx(i)=saxx(i) appy(i)=sayy(i) appz(i)=sazz(i) 6 continue do 7 i=1,3 saxx(i)=appz(i) sayy(i)=appx(i) sazz(i)=appy(i) 7 continue csa=saxx(1)*saxsys(1)+saxx(2)*saxsys(2)+ * saxx(3)*saxsys(3) csb=sayy(1)*saysys(1)+sayy(2)*saysys(2)+ * sayy(3)*saysys(3) csg=sazz(1)*sazsys(1)+sazz(2)*sazsys(2)+ * sazz(3)*sazsys(3) write (*,23) acos(csa)*180./PI,acos(csb)*180./PI,acos(csg)*180./PI write (3,23) acos(csa)*180./PI,acos(csb)*180./PI,acos(csg)*180./PI WRITE(*,'(A)') ' FIVE TENSOR PARAMETERS ' WRITE(3,'(A)') ' ' WRITE(3,'(A)') ' FIVE TENSOR PARAMETERS ' CALL PANG(1,1) b1=-(A1+1.5*A2)/2.0 b2=A1-A2/2.0 WRITE (*,22) * b1,'(',b1/fattoreconv*1e32,' m^3 e+32) ', * b2,'(',b2/fattoreconv*1e32,' m^3 e+32) ' WRITE (3,22) * b1,'(',b1/fattoreconv*1e32,' m^3 e+32) ', * b2,'(',b2/fattoreconv*1e32,' m^3 e+32) ' c c vettori da visualizzare c tsx(1,1)=rccptx(1,2)+saxx(1) tsx(1,2)=rccpty(1,2)+saxx(2) tsx(1,3)=rccptz(1,2)+saxx(3) tsy(1,1)=rccptx(1,2)+sayy(1) tsy(1,2)=rccpty(1,2)+sayy(2) tsy(1,3)=rccptz(1,2)+sayy(3) tsz(1,1)=rccptx(1,2)+sazz(1) tsz(1,2)=rccpty(1,2)+sazz(2) tsz(1,3)=rccptz(1,2)+sazz(3) c write (3,24) saxx(1),saxx(2),saxx(3) c write (3,25) sayy(1),sayy(2),sayy(3) c write (3,26) sazz(1),sazz(2),sazz(3) numberatom=numberatom+1 namatom=' TEX' WRITE(3,'(A)') ' ' WRITE(3,'(A)') ' PDB COORDINATES OF TENSOR SYSTEM' WRITE(3,'(A4,1x,i6,1x,a4,1x,a3,2x,i4,4x,3f8.3)') * 'ATOM',numberatom,namatom,namress,nores, * tsx(1,1),tsx(1,2),tsx(1,3) numberatom=numberatom+1 namatom=' TEY' WRITE(3,'(A4,1x,i6,1x,a4,1x,a3,2x,i4,4x,3f8.3)') * 'ATOM',numberatom,namatom,namress,nores, * tsy(1,1),tsy(1,2),tsy(1,3) numberatom=numberatom+1 namatom=' TEZ' WRITE(3,'(A4,1x,i6,1x,a4,1x,a3,2x,i4,4x,3f8.3)') * 'ATOM',numberatom,namatom,namress,nores, * tsz(1,1),tsz(1,2),tsz(1,3) numberatom=numberatom-3 WRITE(*,'(A,A)') * ' ANGLES BETWEEN REFERENCE AXIS AND TENSOR AXIS ' , * '(degrees) 2 PERMUTATION' WRITE(3,'(A)') ' ' WRITE(3,'(A,A)') * ' ANGLES BETWEEN REFERENCE AXIS AND TENSOR AXIS ', * '(degrees) 2 PERMUTATION' do 8 i=1,3 appx(i)=saxx(i) appy(i)=sayy(i) appz(i)=sazz(i) 8 continue do 9 i=1,3 saxx(i)=appz(i) sayy(i)=appx(i) sazz(i)=appy(i) 9 continue csa=saxx(1)*saxsys(1)+saxx(2)*saxsys(2)+ * saxx(3)*saxsys(3) csb=sayy(1)*saysys(1)+sayy(2)*saysys(2)+ * sayy(3)*saysys(3) csg=sazz(1)*sazsys(1)+sazz(2)*sazsys(2)+ * sazz(3)*sazsys(3) write (*,23) acos(csa)*180./PI,acos(csb)*180./PI,acos(csg)*180./PI write (3,23) acos(csa)*180./PI,acos(csb)*180./PI,acos(csg)*180./PI WRITE(*,'(A)') ' FIVE TENSOR PARAMETERS ' WRITE(3,'(A)') ' ' WRITE(3,'(A)') ' FIVE TENSOR PARAMETERS ' CALL PANG(2,1) b1=-(A1-1.5*A2)/2.0 b2=-(A1+A2/2.0) WRITE (*,22) * b1,'(',b1/fattoreconv*1e32,' m^3 e+32) ', * b2,'(',b2/fattoreconv*1e32,' m^3 e+32) ' WRITE (3,22) * b1,'(',b1/fattoreconv*1e32,' m^3 e+32) ', * b2,'(',b2/fattoreconv*1e32,' m^3 e+32) ' c c vettori da visualizzare c tsx(1,1)=rccptx(1,2)+saxx(1) tsx(1,2)=rccpty(1,2)+saxx(2) tsx(1,3)=rccptz(1,2)+saxx(3) tsy(1,1)=rccptx(1,2)+sayy(1) tsy(1,2)=rccpty(1,2)+sayy(2) tsy(1,3)=rccptz(1,2)+sayy(3) tsz(1,1)=rccptx(1,2)+sazz(1) tsz(1,2)=rccpty(1,2)+sazz(2) tsz(1,3)=rccptz(1,2)+sazz(3) numberatom=numberatom+1 namatom=' TEX' WRITE(3,'(A)') ' ' WRITE(3,'(A)') ' PDB COORDINATES OF TENSOR SYSTEM' WRITE(3,'(A4,1x,i6,1x,a4,1x,a3,2x,i4,4x,3f8.3)') * 'ATOM',numberatom,namatom,namress,nores, * tsx(1,1),tsx(1,2),tsx(1,3) numberatom=numberatom+1 namatom=' TEY' WRITE(3,'(A4,1x,i6,1x,a4,1x,a3,2x,i4,4x,3f8.3)') * 'ATOM',numberatom,namatom,namress,nores, * tsy(1,1),tsy(1,2),tsy(1,3) numberatom=numberatom+1 namatom=' TEZ' WRITE(3,'(A4,1x,i6,1x,a4,1x,a3,2x,i4,4x,3f8.3)') * 'ATOM',numberatom,namatom,namress,nores, * tsz(1,1),tsz(1,2),tsz(1,3) WRITE(3,'(A)') ' ' saxsys(1)=rccptx(1,2)+saxsys(1) saxsys(2)=rccpty(1,2)+saxsys(2) saxsys(3)=rccptz(1,2)+saxsys(3) saysys(1)=rccptx(1,2)+saysys(1) saysys(2)=rccpty(1,2)+saysys(2) saysys(3)=rccptz(1,2)+saysys(3) sazsys(1)=rccptx(1,2)+sazsys(1) sazsys(2)=rccpty(1,2)+sazsys(2) sazsys(3)=rccptz(1,2)+sazsys(3) numberatom=numberatom+1 namatom=' ISX' WRITE(3,'(A)') ' ' WRITE(3,'(A)') ' PDB COORDINATES OF REFERENCE SYSTEM' WRITE(3,'(A4,1x,i6,1x,a4,1x,a3,2x,i4,4x,3f8.3)') * 'ATOM',numberatom,namatom,namress,nores, * saxsys(1),saxsys(2),saxsys(3) numberatom=numberatom+1 namatom=' ISY' WRITE(3,'(A4,1x,i6,1x,a4,1x,a3,2x,i4,4x,3f8.3)') * 'ATOM',numberatom,namatom,namress,nores, * saysys(1),saysys(2),saysys(3) numberatom=numberatom+1 namatom=' ISZ' WRITE(3,'(A4,1x,i6,1x,a4,1x,a3,2x,i4,4x,3f8.3)') * 'ATOM',numberatom,namatom,namress,nores, * sazsys(1),sazsys(2),sazsys(3) 10 continue close(2) 22 FORMAT * (1X,'A1= ',F9.3,2x,a1,F9.4,a,2X,'A2= ',F9.3,2x,a1,F9.4,a) 23 FORMAT (1X,'arccos x^xt = ',F8.3,4X,'arccos y^yt = ',F8.3, * 4X,'arccos z^zt = ',F8.3) c 24 format (1X,'assex : ',1X,F12.3,1X,F12.3,1X,F12.3) c 25 format (1X,'assey : ',1X,F12.3,1X,F12.3,1X,F12.3) c 26 format (1X,'assez : ',1X,F12.3,1X,F12.3,1X,F12.3) END C-------------------------------------------------------------- C C PANG C C This subroutine computes the three new angles of Euler for C each permutation of tensor axes C-------------------------------------------------------------- SUBROUTINE PANG(n,nr) C c n = Number of axes permutation: c 0 = none; 1 = first permutation; else = second perm. c nr = Number (label) of reference system c include 'fantacross.h' dimension temp(3) ! Auxiliary variable for permutations PI=3.1415927 c Composing permuted rotation matrix: c Note: asxx(i,j) = 1st column of rot. matrix of tensor i c asyy(i,j) = 2nd column of rot. matrix of tensor i c aszz(i,j) = 3rd column of rot. matrix of tensor i c (j=1,2,3; row index) c----------------------------------------------- if (n.eq.0) then ! case 0: No permutation c x is x sxx=asxx(nr,1) sxy=asxx(nr,2) sxz=asxx(nr,3) c y is y syx=asyy(nr,1) syy=asyy(nr,2) syz=asyy(nr,3) c z is z szx=aszz(nr,1) szy=aszz(nr,2) szz=aszz(nr,3) elseif (n.eq.1) then ! case 1: First permutation of axes c x is z sxx=aszz(nr,1) sxy=aszz(nr,2) sxz=aszz(nr,3) c y is x syx=asxx(nr,1) syy=asxx(nr,2) syz=asxx(nr,3) c z is y szx=asyy(nr,1) szy=asyy(nr,2) szz=asyy(nr,3) else ! case default: Second permutation of axes c x is y sxx=asyy(nr,1) sxy=asyy(nr,2) sxz=asyy(nr,3) c y is z syx=aszz(nr,1) syy=aszz(nr,2) syz=aszz(nr,3) c z is x szx=asxx(nr,1) szy=asxx(nr,2) szz=asxx(nr,3) endif c Computing angles of Euler: c let i,j,k = versors of global reference system; c sx,sy,sz = versors of tensor; c and N = (k) x (sx) the Node line; then: c c phi is the angle from j to N (right-hand rotation on k) c theta is the angle from N to sy (right-hand rotation on sx) c omega is the angle from ij plane to sx (left-hand rot. on N). c c In case of sx=k then N is null and: c c phi+theta is the angle from j to sy (right-hand rotation on k) c omega is +90 (deg) c (the single values of phi and theta doesn't matter, so we can c assume theta=0 and phi=phi+theta). c c In case of sx=-k then N is null and: c phi-theta is the angle from j to sy (right-hand rotation on k) c omega is -90 (deg) c (the single values of phi and theta doesn't matter, so we can c assume theta=0 and phi=phi-theta). c c----------------------------------------------- c Node line: N=(k)x(sx)= -sxy i + sxx j absN=sqrt(sxy**2+sxx**2) ! Node line module if(absN.lt.EPS) then ! N is null and sx is parallel to k t=0 if(sxz.gt.0) then ! sx=k o=PI/2. else ! sx=-k o=-PI/2. endif p=acos(syy) ! cos(p)=(sy)(j) [scalar product] = syy c Phi sign correction: c sign must switch if (j x sy)(k)<0 if(syx.gt.0) p=-p ! rotation on k was left-hand else c Phi: p=acos(sxx/absN) c sign must switch if (j x N)(k)<0 if(sxy.lt.0) p=-p ! rotation on k was left-hand c Theta: t=acos((syy*sxx-syx*sxy)/absN) c sign must switch if (N x sy)(sz)<0 c computing (N x sy)(sz): scalpr=(sxx*sxx+sxy*sxy)*syz-sxz*(sxy*syy+sxx*syx) if(scalpr.lt.0) t=-t ! rotation on sx was left-hand c Omega: o=acos(absN) c sign must switch if (sx x sx')(N)<0 [where sx' is the c projection of sx on plane xy], that is: c if (sxx**2+sxy**2)*sxz > 0 ---> sxz > 0 if(sxz.lt.0) o=-o ! rotation on N was right-hand endif WRITE (*,21) p,t,o WRITE (3,21) p,t,o 21 FORMAT(1X,'PHI= ',F7.4,4X,'THETA= ',F7.4,4X,'OMEGA= ',F7.4) END SUBROUTINE F01BCF(N,TOL,Z,IZ,W,IW,D,E,C,S) C MARK 3 RELEASE. NAG COPYRIGHT 1972. C MARK 4 REVISED. C MARK 4.5 REVISED C MARK 5C REVISED C MARK 6B REVISED IER-125 IER-127 (JUL 1978) C MARK 11 REVISED. VECTORISATION (JAN 1984). C MARK 11.5(F77) REVISED. (SEPT 1985.) C C C TRECX2 C F01BCF REDUCES A COMPLEX HERMITIAN MATRIX TO REAL C TRIDIAGONAL FORM FROM WHICH THE EIGENVALUES AND EIGENVECTORS C CAN BE FOUND USING SUBROUTINE F02AYF,(CXTQL2). THE HERMITIAN C MATRIX A=A(1) IS REDUCED TO THE TRIDIAGONAL MATRIX A(N-1) BY C N-2 UNITARY TRANSFORMATIONS. THE HOUSEHOLDER REDUCTION ITSELF C DOES NOT GIVE A REAL TRIDIAGONAL MATRIX, THE OFF-DIAGONAL C ELEMENTS ARE COMPLEX. THEY ARE SUBSEQUENTLY MADE REAL BY A C DIAGONAL TRANSFORMATION. C APRIL 1ST. 1972 C C .. SCALAR ARGUMENTS .. DOUBLE PRECISION TOL INTEGER IW, IZ, N C .. ARRAY ARGUMENTS .. DOUBLE PRECISION C(N), D(N), E(N), S(N), W(IW,N), Z(IZ,N) C .. LOCAL SCALARS .. DOUBLE PRECISION CO, F, FI, FR, G, GI, GR, H, HH, R, SI INTEGER I, II, J, K, L C .. EXTERNAL SUBROUTINES .. EXTERNAL F01BCY, F01BCZ C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS, SQRT C .. EXECUTABLE STATEMENTS .. DO 20 I = 1, N D(I) = Z(N,I) E(I) = -W(N,I) 20 CONTINUE IF (N.EQ.1) GO TO 540 DO 360 II = 2, N I = N - II + 2 L = I - 2 G = 0.0D0 FR = D(I-1) FI = E(I-1) IF (L.EQ.0) GO TO 60 DO 40 K = 1, L G = G + D(K)*D(K) + E(K)*E(K) 40 CONTINUE 60 H = G + FR*FR + FI*FI C L IS NOW I-1 L = L + 1 IF (ABS(FR)+ABS(FI).NE.0.0D0) GO TO 80 R = 0.0D0 CO = 1.0D0 C(I) = 1.0D0 SI = 0.0D0 S(I) = 0.0D0 GO TO 140 80 IF (ABS(FR).LT.ABS(FI)) GO TO 100 R = ABS(FR)*SQRT(1.0D0+(FI/FR)**2) GO TO 120 100 R = ABS(FI)*SQRT(1.0D0+(FR/FI)**2) 120 SI = FI/R S(I) = -SI CO = FR/R C(I) = CO 140 IF (G.LE.TOL) GO TO 280 G = -SQRT(H) E(I) = G C E(I) HAS ITS FINAL REAL VALUE H = H - R*G C S*S + SR D(I-1) = (R-G)*CO E(I-1) = (R-G)*SI DO 160 J = 1, L Z(J,I) = D(J) W(J,I) = E(J) 160 CONTINUE CALL F01BCZ(Z,IZ,W,IW,L,D,E,C,S) C FORM P DO 180 J = 1, L C(J) = C(J)/H S(J) = S(J)/H 180 CONTINUE FR = 0.0D0 DO 200 J = 1, L FR = FR + C(J)*D(J) + S(J)*E(J) 200 CONTINUE C FORM K HH = FR/(H+H) C FORM Q DO 220 J = 1, L C(J) = C(J) - HH*D(J) S(J) = S(J) - HH*E(J) 220 CONTINUE C NOW FORM REDUCED A DO 260 J = 1, L FR = D(J) FI = E(J) GR = C(J) GI = S(J) DO 240 K = J, L Z(K,J) = (((Z(K,J)-GR*D(K))-GI*E(K))-FR*C(K)) - FI*S(K) W(K,J) = (((W(K,J)-GR*E(K))+GI*D(K))-FR*S(K)) + FI*C(K) 240 CONTINUE D(J) = Z(L,J) Z(I,J) = 0.0D0 E(J) = -W(L,J) W(I,J) = 0.0D0 W(J,J) = 0.0D0 260 CONTINUE GO TO 340 280 E(I) = R H = 0.0D0 DO 300 J = 1, L Z(J,I) = D(J) W(J,I) = E(J) 300 CONTINUE DO 320 J = 1, L Z(I,J) = 0.0D0 D(J) = Z(I-1,J) W(I,J) = 0.0D0 E(J) = -W(I-1,J) 320 CONTINUE 340 D(I) = H 360 CONTINUE C WE NOW FORM THE PRODUCT OF THE C HOUSEHOLDER MATRICES, OVERWRITING C ON Z AND W DO 500 I = 2, N L = I - 1 Z(N,L) = Z(L,L) Z(L,L) = 1.0D0 W(N,L) = E(L) W(L,L) = 0.0D0 H = D(I) IF (H.EQ.0.0D0) GO TO 460 DO 380 K = 1, L D(K) = 0.0D0 E(K) = 0.0D0 380 CONTINUE CALL F01BCY(Z,IZ,W,IW,L,L,Z(1,I),W(1,I),D,E) DO 400 K = 1, L D(K) = D(K)/H E(K) = -E(K)/H 400 CONTINUE DO 440 J = 1, L DO 420 K = 1, L Z(K,J) = Z(K,J) - Z(K,I)*D(J) + W(K,I)*E(J) W(K,J) = W(K,J) - Z(K,I)*E(J) - W(K,I)*D(J) 420 CONTINUE 440 CONTINUE 460 DO 480 J = 1, L Z(J,I) = 0.0D0 W(J,I) = 0.0D0 480 CONTINUE 500 CONTINUE W(N,N) = E(N) DO 520 I = 1, N D(I) = Z(N,I) Z(N,I) = 0.0D0 E(I) = W(N,I) W(N,I) = 0.0D0 520 CONTINUE 540 Z(N,N) = 1.0D0 W(N,N) = 0.0D0 E(1) = 0.0D0 C NOW WE MULTIPLY BY THE C COSTHETA + I SINTHETA COLUMN C FACTORS CO = 1.0D0 SI = 0.0D0 IF (N.EQ.1) RETURN DO 580 I = 2, N F = CO*C(I) - SI*S(I) SI = CO*S(I) + SI*C(I) CO = F DO 560 J = 1, N F = Z(J,I)*CO - W(J,I)*SI W(J,I) = Z(J,I)*SI + W(J,I)*CO Z(J,I) = F 560 CONTINUE 580 CONTINUE RETURN END SUBROUTINE F01BCY(AR,IAR,AI,IAI,M,N,BR,BI,CR,CI) C MARK 11 RELEASE. NAG COPYRIGHT 1983. C MARK 11.5(F77) REVISED. (SEPT 1985.) C C COMPUTES C = C + (A**H)*B (COMPLEX) WHERE C A IS RECTANGULAR M BY N. C C MUST BE DISTINCT FROM B. C C C .. SCALAR ARGUMENTS .. INTEGER IAI, IAR, M, N C .. ARRAY ARGUMENTS .. DOUBLE PRECISION AI(IAI,N), AR(IAR,N), BI(M), BR(M), CI(N), CR(N) C .. LOCAL SCALARS .. DOUBLE PRECISION XI, XR INTEGER I, J C .. EXECUTABLE STATEMENTS .. DO 40 I = 1, N XR = CR(I) XI = CI(I) DO 20 J = 1, M XR = XR + AR(J,I)*BR(J) + AI(J,I)*BI(J) XI = XI + AR(J,I)*BI(J) - AI(J,I)*BR(J) 20 CONTINUE CR(I) = XR CI(I) = XI 40 CONTINUE RETURN END SUBROUTINE F01BCZ(AR,IAR,AI,IAI,N,BR,BI,CR,CI) C MARK 11 RELEASE. NAG COPYRIGHT 1983. C MARK 11.5(F77) REVISED. (SEPT 1985.) C C COMPUTES C = A*B (COMPLEX) WHERE C A IS A HERMITIAN N-BY-N MATRIX, C WHOSE LOWER TRIANGLE IS STORED IN A. C C MUST BE DISTINCT FROM B. C C C .. SCALAR ARGUMENTS .. INTEGER IAI, IAR, N C .. ARRAY ARGUMENTS .. DOUBLE PRECISION AI(IAI,N), AR(IAR,N), BI(N), BR(N), CI(N), CR(N) C .. LOCAL SCALARS .. DOUBLE PRECISION YI, YR INTEGER I, IP1, J, NM1 C .. EXECUTABLE STATEMENTS .. DO 20 I = 1, N CR(I) = 0.0D0 CI(I) = 0.0D0 20 CONTINUE IF (N.EQ.1) GO TO 100 NM1 = N - 1 DO 80 I = 1, NM1 DO 40 J = I, N CR(J) = CR(J) + AR(J,I)*BR(I) - AI(J,I)*BI(I) CI(J) = CI(J) + AR(J,I)*BI(I) + AI(J,I)*BR(I) 40 CONTINUE YR = CR(I) YI = CI(I) IP1 = I + 1 DO 60 J = IP1, N YR = YR + AR(J,I)*BR(J) + AI(J,I)*BI(J) YI = YI + AR(J,I)*BI(J) - AI(J,I)*BR(J) 60 CONTINUE CR(I) = YR CI(I) = YI 80 CONTINUE 100 CR(N) = CR(N) + AR(N,N)*BR(N) - AI(N,N)*BI(N) CI(N) = CI(N) + AR(N,N)*BI(N) + AI(N,N)*BR(N) RETURN END SUBROUTINE F02AXF(AR,IAR,AI,IAI,N,WR,VR,IVR,VI,IVI,WK1,WK2,WK3, * IFAIL) C MARK 3 RELEASE. NAG COPYRIGHT 1972. C MARK 4.5 REVISED C MARK 9 REVISED. IER-327 (SEP 1981). C MARK 11.5(F77) REVISED. (SEPT 1985.) C MARK 13 REVISED. USE OF MARK 12 X02 FUNCTIONS (APR 1988). C C EIGENVALUES AND EIGENVECTORS OF A COMPLEX HERMITIAN MATRIX C 1ST APRIL 1972 C C .. PARAMETERS .. CHARACTER*6 SRNAME PARAMETER (SRNAME='F02AXF') C .. SCALAR ARGUMENTS .. INTEGER IAI, IAR, IFAIL, IVI, IVR, N C .. ARRAY ARGUMENTS .. DOUBLE PRECISION AI(IAI,N), AR(IAR,N), VI(IVI,N), VR(IVR,N), * WK1(N), WK2(N), WK3(N), WR(N) C .. LOCAL SCALARS .. DOUBLE PRECISION A, B, MAX, SQ, SUM, XXXX INTEGER I, ISAVE, J C .. LOCAL ARRAYS .. CHARACTER*1 P01REC(1) C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION X02AJF, X02AKF INTEGER P01ABF EXTERNAL X02AJF, X02AKF, P01ABF C .. EXTERNAL SUBROUTINES .. EXTERNAL F01BCF, F02AYF C .. INTRINSIC FUNCTIONS .. INTRINSIC SQRT C .. EXECUTABLE STATEMENTS .. ISAVE = IFAIL DO 40 I = 1, N IF (AI(I,I).NE.0.0D0) GO TO 140 DO 20 J = 1, I VR(I,J) = AR(I,J) VI(I,J) = AI(I,J) 20 CONTINUE 40 CONTINUE CALL F01BCF(N,X02AKF()/X02AJF(),VR,IVR,VI,IVI,WR,WK1,WK2,WK3) IFAIL = 1 CALL F02AYF(N,X02AJF(),WR,WK1,VR,IVR,VI,IVI,IFAIL) IF (IFAIL.EQ.0) GO TO 60 IFAIL = P01ABF(ISAVE,1,SRNAME,0,P01REC) RETURN C NORMALISE 60 DO 120 I = 1, N SUM = 0.0D0 MAX = 0.0D0 DO 80 J = 1, N SQ = VR(J,I)*VR(J,I) + VI(J,I)*VI(J,I) SUM = SUM + SQ IF (SQ.LE.MAX) GO TO 80 MAX = SQ A = VR(J,I) B = VI(J,I) 80 CONTINUE IF (SUM.EQ.0.0D0) GO TO 120 SUM = 1.0D0/SQRT(SUM*MAX) DO 100 J = 1, N SQ = SUM*(VR(J,I)*A+VI(J,I)*B) VI(J,I) = SUM*(VI(J,I)*A-VR(J,I)*B) VR(J,I) = SQ 100 CONTINUE 120 CONTINUE RETURN 140 IFAIL = P01ABF(ISAVE,2,SRNAME,0,P01REC) RETURN END SUBROUTINE F02AYF(N,EPS,D,E,Z,IZ,W,IW,IFAIL) C MARK 3 RELEASE. NAG COPYRIGHT 1972. C MARK 4 REVISED. C MARK 4.5 REVISED C MARK 9 REVISED. IER-326 (SEP 1981). C MARK 11.5(F77) REVISED. (SEPT 1985.) C C CXTQL2 C THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS OF A C HERMITIAN MATRIX, WHICH HAS BEEN REDUCED TO A REAL C TRIDIAGONAL MATRIX, T, GIVEN WITH ITS DIAGONAL ELEMENTS IN C THE ARRAY D(N) AND ITS SUB-DIAGONAL ELEMENTS IN THE LAST N C - 1 STORES OF THE ARRAY E(N), USING QL TRANSFORMATIONS. THE C EIGENVALUES ARE OVERWRITTEN ON THE DIAGONAL ELEMENTS IN THE C ARRAY D IN ASCENDING ORDER. THE REAL AND IMAGINARY PARTS OF C THE EIGENVECTORS ARE FORMED IN THE ARRAYS Z,W(N,N) C RESPECTIVELY, OVERWRITING THE ACCUMULATED TRANSFORMATIONS AS C SUPPLIED BY THE SUBROUTINE F01BCF. THE SUBROUTINE WILL FAIL C IF ALL EIGENVALUES TAKE MORE THAN 30*N ITERATIONS C 1ST APRIL 1972 C C .. PARAMETERS .. CHARACTER*6 SRNAME PARAMETER (SRNAME='F02AYF') C .. SCALAR ARGUMENTS .. DOUBLE PRECISION EPS INTEGER IFAIL, IW, IZ, N C .. ARRAY ARGUMENTS .. DOUBLE PRECISION D(N), E(N), W(IW,N), Z(IZ,N) C .. LOCAL SCALARS .. DOUBLE PRECISION B, C, F, G, H, P, R, S INTEGER I, I1, II, ISAVE, J, K, L, M, M1 C .. LOCAL ARRAYS .. CHARACTER*1 P01REC(1) C .. EXTERNAL FUNCTIONS .. INTEGER P01ABF EXTERNAL P01ABF C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS, SQRT C .. EXECUTABLE STATEMENTS .. ISAVE = IFAIL IF (N.EQ.1) GO TO 40 DO 20 I = 2, N E(I-1) = E(I) 20 CONTINUE 40 E(N) = 0.0D0 B = 0.0D0 F = 0.0D0 J = 30*N DO 300 L = 1, N H = EPS*(ABS(D(L))+ABS(E(L))) IF (B.LT.H) B = H C LOOK FOR SMALL SUB-DIAG ELEMENT DO 60 M = L, N IF (ABS(E(M)).LE.B) GO TO 80 60 CONTINUE 80 IF (M.EQ.L) GO TO 280 100 IF (J.LE.0) GO TO 400 J = J - 1 C FORM SHIFT G = D(L) H = D(L+1) - G IF (ABS(H).GE.ABS(E(L))) GO TO 120 P = H*0.5D0/E(L) R = SQRT(P*P+1.0D0) H = P + R IF (P.LT.0.0D0) H = P - R D(L) = E(L)/H GO TO 140 120 P = 2.0D0*E(L)/H R = SQRT(P*P+1.0D0) D(L) = E(L)*P/(1.0D0+R) 140 H = G - D(L) I1 = L + 1 IF (I1.GT.N) GO TO 180 DO 160 I = I1, N D(I) = D(I) - H 160 CONTINUE 180 F = F + H C QL TRANSFORMATION P = D(M) C = 1.0D0 S = 0.0D0 M1 = M - 1 DO 260 II = L, M1 I = M1 - II + L G = C*E(I) H = C*P IF (ABS(P).LT.ABS(E(I))) GO TO 200 C = E(I)/P R = SQRT(C*C+1.0D0) E(I+1) = S*P*R S = C/R C = 1.0D0/R GO TO 220 200 C = P/E(I) R = SQRT(C*C+1.0D0) E(I+1) = S*E(I)*R S = 1.0D0/R C = C/R 220 P = C*D(I) - S*G D(I+1) = H + S*(C*G+S*D(I)) C FORM VECTOR DO 240 K = 1, N H = Z(K,I+1) Z(K,I+1) = S*Z(K,I) + C*H Z(K,I) = C*Z(K,I) - S*H H = W(K,I+1) W(K,I+1) = S*W(K,I) + C*H W(K,I) = C*W(K,I) - S*H 240 CONTINUE 260 CONTINUE E(L) = S*P D(L) = C*P IF (ABS(E(L)).GT.B) GO TO 100 280 D(L) = D(L) + F 300 CONTINUE C ORDER EIGEN VALUES ANS EIGENVECTORS DO 380 I = 1, N K = I P = D(I) I1 = I + 1 IF (I1.GT.N) GO TO 340 DO 320 J = I1, N IF (D(J).GE.P) GO TO 320 K = J P = D(J) 320 CONTINUE 340 IF (K.EQ.I) GO TO 380 D(K) = D(I) D(I) = P DO 360 J = 1, N P = Z(J,I) Z(J,I) = Z(J,K) Z(J,K) = P P = W(J,I) W(J,I) = W(J,K) W(J,K) = P 360 CONTINUE 380 CONTINUE IFAIL = 0 RETURN 400 IFAIL = P01ABF(ISAVE,1,SRNAME,0,P01REC) RETURN END INTEGER FUNCTION P01ABF(IFAIL,IERROR,SRNAME,NREC,REC) C MARK 11.5(F77) RELEASE. NAG COPYRIGHT 1986. C MARK 13 REVISED. IER-621 (APR 1988). C MARK 13B REVISED. IER-668 (AUG 1988). C C P01ABF IS THE ERROR-HANDLING ROUTINE FOR THE NAG LIBRARY. C C P01ABF EITHER RETURNS THE VALUE OF IERROR THROUGH THE ROUTINE C NAME (SOFT FAILURE), OR TERMINATES EXECUTION OF THE PROGRAM C (HARD FAILURE). DIAGNOSTIC MESSAGES MAY BE OUTPUT. C C IF IERROR = 0 (SUCCESSFUL EXIT FROM THE CALLING ROUTINE), C THE VALUE 0 IS RETURNED THROUGH THE ROUTINE NAME, AND NO C MESSAGE IS OUTPUT C C IF IERROR IS NON-ZERO (ABNORMAL EXIT FROM THE CALLING ROUTINE), C THE ACTION TAKEN DEPENDS ON THE VALUE OF IFAIL. C C IFAIL = 1: SOFT FAILURE, SILENT EXIT (I.E. NO MESSAGES ARE C OUTPUT) C IFAIL = -1: SOFT FAILURE, NOISY EXIT (I.E. MESSAGES ARE OUTPUT) C IFAIL =-13: SOFT FAILURE, NOISY EXIT BUT STANDARD MESSAGES FROM C P01ABF ARE SUPPRESSED C IFAIL = 0: HARD FAILURE, NOISY EXIT C C FOR COMPATIBILITY WITH CERTAIN ROUTINES INCLUDED BEFORE MARK 12 C P01ABF ALSO ALLOWS AN ALTERNATIVE SPECIFICATION OF IFAIL IN WHICH C IT IS REGARDED AS A DECIMAL INTEGER WITH LEAST SIGNIFICANT DIGITS C CBA. THEN C C A = 0: HARD FAILURE A = 1: SOFT FAILURE C B = 0: SILENT EXIT B = 1: NOISY EXIT C C EXCEPT THAT HARD FAILURE NOW ALWAYS IMPLIES A NOISY EXIT. C C S.HAMMARLING, M.P.HOOPER AND J.J.DU CROZ, NAG CENTRAL OFFICE. C C .. SCALAR ARGUMENTS .. INTEGER IERROR, IFAIL, NREC CHARACTER*(*) SRNAME C .. ARRAY ARGUMENTS .. CHARACTER*(*) REC(*) C .. LOCAL SCALARS .. INTEGER I, NERR CHARACTER*72 MESS C .. EXTERNAL SUBROUTINES .. EXTERNAL P01ABZ, X04AAF, X04BAF C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS, MOD C .. EXECUTABLE STATEMENTS .. IF (IERROR.NE.0) THEN C ABNORMAL EXIT FROM CALLING ROUTINE IF (IFAIL.EQ.-1 .OR. IFAIL.EQ.0 .OR. IFAIL.EQ.-13 .OR. * (IFAIL.GT.0 .AND. MOD(IFAIL/10,10).NE.0)) THEN C NOISY EXIT CALL X04AAF(0,NERR) DO 20 I = 1, NREC CALL X04BAF(NERR,REC(I)) 20 CONTINUE IF (IFAIL.NE.-13) THEN WRITE (MESS,FMT=99999) SRNAME, IERROR CALL X04BAF(NERR,MESS) IF (ABS(MOD(IFAIL,10)).NE.1) THEN C HARD FAILURE CALL X04BAF(NERR, * ' ** NAG HARD FAILURE - EXECUTION TERMINATED' * ) CALL P01ABZ ELSE C SOFT FAILURE CALL X04BAF(NERR, * ' ** NAG SOFT FAILURE - CONTROL RETURNED') END IF END IF END IF END IF P01ABF = IERROR RETURN C 99999 FORMAT (' ** ABNORMAL EXIT FROM NAG LIBRARY ROUTINE ',A,': IFAIL', * ' =',I6) END SUBROUTINE P01ABZ C MARK 11.5(F77) RELEASE. NAG COPYRIGHT 1986. C C TERMINATES EXECUTION WHEN A HARD FAILURE OCCURS. C C ******************** IMPLEMENTATION NOTE ******************** C THE FOLLOWING STOP STATEMENT MAY BE REPLACED BY A CALL TO AN C IMPLEMENTATION-DEPENDENT ROUTINE TO DISPLAY A MESSAGE AND/OR C TO ABORT THE PROGRAM. C ************************************************************* C .. EXECUTABLE STATEMENTS .. STOP END DOUBLE PRECISION FUNCTION X02AJF() C MARK 12 RELEASE. NAG COPYRIGHT 1986. C C RETURNS (1/2)*B**(1-P) IF ROUNDS IS .TRUE. C RETURNS B**(1-P) OTHERWISE C C .. EXECUTABLE STATEMENTS .. X02AJF = 0.11102230246251568000D-015 RETURN END DOUBLE PRECISION FUNCTION X02AKF() C MARK 12 RELEASE. NAG COPYRIGHT 1986. C C RETURNS B**(EMIN-1) (THE SMALLEST POSITIVE MODEL NUMBER) C C .. EXECUTABLE STATEMENTS .. X02AKF = 0.22250738585072014000D-307 RETURN END SUBROUTINE X04AAF(I,NERR) C MARK 7 RELEASE. NAG COPYRIGHT 1978 C MARK 7C REVISED IER-190 (MAY 1979) C MARK 11.5(F77) REVISED. (SEPT 1985.) C IF I = 0, SETS NERR TO CURRENT ERROR MESSAGE UNIT NUMBER C (STORED IN NERR1). C IF I = 1, CHANGES CURRENT ERROR MESSAGE UNIT NUMBER TO C VALUE SPECIFIED BY NERR. C C .. SCALAR ARGUMENTS .. INTEGER I, NERR C .. LOCAL SCALARS .. INTEGER NERR1 C .. SAVE STATEMENT .. SAVE NERR1 C .. DATA STATEMENTS .. DATA NERR1/0/ C .. EXECUTABLE STATEMENTS .. IF (I.EQ.0) NERR = NERR1 IF (I.EQ.1) NERR1 = NERR RETURN END SUBROUTINE X04BAF(NOUT,REC) C MARK 11.5(F77) RELEASE. NAG COPYRIGHT 1986. C C X04BAF WRITES THE CONTENTS OF REC TO THE UNIT DEFINED BY NOUT. C C TRAILING BLANKS ARE NOT OUTPUT, EXCEPT THAT IF REC IS ENTIRELY C BLANK, A SINGLE BLANK CHARACTER IS OUTPUT. C IF NOUT.LT.0, I.E. IF NOUT IS NOT A VALID FORTRAN UNIT IDENTIFIER, C THEN NO OUTPUT OCCURS. C C .. SCALAR ARGUMENTS .. INTEGER NOUT CHARACTER*(*) REC C .. LOCAL SCALARS .. INTEGER I C .. INTRINSIC FUNCTIONS .. INTRINSIC LEN C .. EXECUTABLE STATEMENTS .. IF (NOUT.GE.0) THEN C REMOVE TRAILING BLANKS DO 20 I = LEN(REC), 2, -1 IF (REC(I:I).NE.' ') GO TO 40 20 CONTINUE C WRITE RECORD TO EXTERNAL FILE 40 WRITE (NOUT,FMT=99999) REC(1:I) END IF RETURN C 99999 FORMAT (A) END