PROGRAM HYDR IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (MX3=3*1100,MX=MX3/3,n40=4) PARAMETER (NG0=100) DIMENSION RVAC(MX3),RSOL(MX3),ITVAC(MX3),ITSOL(MX3), 1FVAC(MX3,MX3),IND(MX),pot(n40),p(n40),e(3),IBT(7*MX), 2iagv(NG0),jagv(NG0),kagv(NG0),lagv(NG0), 3iags(NG0),jags(NG0),kags(NG0),lags(NG0),il(MX),im(MX), 4PA(MX,3,3),A(MX,3,3),AJ(MX,3,3),PV(MX,3,3),pd(n40) data p/-0.00325d0,-0.00624d0,0.00269d0,0.00680/ data pd/0.970d0,-1.842d0,0.209d0,0.663d0/ data w0/1731.72d0/ data d0/0.0511199d0/ LOGICAL lcct,lten,lper C unitf=1.602D-19**2/(4.0d0*3.141592d0*8.854D-12)/ 1 1.0D-10*2.1066D26/4.186D3 c OPEN(16,FILE='FTRY.OUT') WRITE(*,60000) WRITE(16,60000) 60000 FORMAT(' EMPIRICAL SOLVENT (WATER) FORCE FIELD CORRECTION',/,/, 3 ' Petr Bour, Prague/Chicago 2003',/,/, 1 ' Input files : VAC.X (SMALL)' ,/, 1 ' VAC.FC',/, 1 ' SOL.XYZ (BIG)',/, 1 ' CCT.INP (OPTIONAL, OVERLAP)',/, 1 ' Output files: FILE.X' ,/, 1 ' FILE.FC',/,/) C c INQUIRE(FILE='PER.KEY',exist=lper) AAXIS=0.0d0 if(lper)then open(15,file='PER.KEY') read(15,*)AAXIS CLOSE(15) write(6,6700)AAXIS 6700 format(' Periodic boundary axis ',F10.3) endif INQUIRE(FILE='VAC.TEN',exist=lten) if(lten)call READTEN(MX,PA,A,AJ,PV,NOPT,NQ) OPEN(35,FILE='VAC.X',STATUS='OLD') CALL READX(RVAC,ITVAC,NATVAC,IBT) WRITE(*,*)' Vacuum geometry read in ... ' OPEN(35,FILE='SOL.XYZ',STATUS='OLD') CALL READXYZ(RSOL,ITSOL,NATSOL) WRITE(*,*)' Solvent geometry read in ... ' OPEN(20,FILE='VAC.FC',STATUS='OLD') CALL READFF(MX3,NATVAC*3,FVAC) WRITE(*,*)' Vacuum FF read in ... ' pisvejc=1.1191D0 pis2=0.72d0 WRITE(6,60777)pisvejc,pis2 60777 format(' Pisvejc',1h','s constants: ',2f8.4) C c Look for C-(C=O)-N group call findag(NATVAC,RVAC,ITVAC,iagv,jagv,kagv,lagv,igv,NG0) call findag(NATSOL,RSOL,ITSOL,iags,jags,kags,lags,igs,NG0) if(igv.ne.igs)then write(6,*)'Different number of AGs!' write(6,*)igv,igs stop endif write(6,799)igv 799 format(i5,' amide groups',/, 1' group vacuum solvent') do 7 i=1,igv 7 write(6,800)i,iagv(i),jagv(i),kagv(i),lagv(i), 1 iags(i),jags(i),kags(i),lags(i) 800 format(i5,4i5,5x,4i5) c call rCCT(NATSOL,IND,lcct) if(lcct)then WRITE(*,*)' Overlap read in ... ' OPEN(34,file='BIG.X') write(34,*) write(34,*)NATVAC do 90 ia=1,NATVAC iroz=0 do 91 ii=1,NATSOL 91 if(IND(ii).eq.ia)iroz=ii if(iroz.eq.0)then write(6,*)'iroz=0!' stop endif 90 write(34,3400)ITVAC(ia),(RSOL(3*(iroz-1)+ix),ix=1,3), 1(IBT(7*(ia-1)+jj),jj=1,7) 3400 format(i5,3f15.7,' ',7i6,' 0.0') close(34) WRITE(*,*)' BIG.X written ...' else WRITE(*,*)' Overlap not found, automatic assignement used' call ass(iags,jags,kags,lags,iagv,jagv,kagv,lagv, 1ind,il,im,RVAC,RSOL,igs,igv) endif do 1 ig=1,igv i=iagv(ig) j=jagv(ig) k=kagv(ig) l=lagv(ig) write(6,6001)ig,k,j,i,l 6001 format(i4,'. group vac C (C = O ) N :',4i10) ii=0 jj=0 kk=0 ll=0 do 5 iind=1,NATSOL if(ind(iind).eq.i)ii=iind if(ind(iind).eq.j)jj=iind if(ind(iind).eq.k)kk=iind 5 if(ind(iind).eq.l)ll=iind if(ii.ne.0.and.jj.ne.0.and.kk.ne.0.and.ll.ne.0)then call getpot(RSOL,ii,jj,kk,ll,pot,NATSOL,ITSOL,AAXIS,lper) delw=0.0d0 delD=0.0d0 do 6 ia=1,4 delD=delD+pot(ia)*pd(ia) 6 delw=delw+pot(ia)* p(ia)*unitf write(6,6002)kk,jj,ii,ll,pot(3),pot(2),pot(1),pot(4), 1 p(3),p(2),p(1),p(4),pd(3),pd(2),pd(1),pd(4),w0+delw, 2d0+delD 6002 format(4x,' sol C (C = O ) N :',4i10,/, 1 4x,' potential :',4f10.5,/, 2 4x,' parameters :',4f10.5,/, 2 4x,' D parameters :',4f10.5,/, 3 4x,' w :',f10.2,/, 3 4x,' d :',f10.5) c c Repair force field: call repair(FVAC,i,j,delw,w0,MX3,RVAC,pisvejc,e) if(lten)call repait(PA,i,j,delw,w0,MX,pis2,e,d0,delD) else write(6,6003) 6003 format(4x,' sol not found') stop endif 1 continue call WRITEFF(MX3,3*NATVAC,FVAC) WRITE(*,*)' FILE.FC written' OPEN(35,FILE='FILE.X') call WRITEX(RVAC,ITVAC,NATVAC) WRITE(*,*)' FILE.X written' if(lten)call WRITETEN(MX,PA,A,AJ,PV,NOPT,NQ,NATVAC) WRITE(*,*)' << Program finished OK >>' STOP END C SUBROUTINE WRITEFF(MX3,N,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(MX3,N) OPEN(20,FILE='FILE.FC') N1=1 1 N3=N1+4 IF(N3.GT.N)N3=N DO 130 LN=N1,N 130 WRITE(20,17)LN,(FCAR(LN,J),J=N1,MIN(LN,N3)) N1=N1+5 IF(N3.LT.N)GOTO 1 17 FORMAT(I4,5D14.6) CLOSE(20) RETURN END c SUBROUTINE READXYZ(R,IT,NAT) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8(A-H,O-Z) DIMENSION R(*),IT(*) character*2 at READ(35,*)NAT DO 3 I=1,NAT READ(35,355)idum,at,(R(3*(I-1)+j),j=1,3) 355 format(i6,2x,a2,1x,3f12.6) IT(I)=0 if(at(1:1).eq.'H')IT(I)=1 if(at(1:1).eq.'C')IT(I)=6 if(at(1:1).eq.'N')IT(I)=7 if(at(1:1).eq.'O')IT(I)=8 if(IT(I).eq.0)write(6,*)'unknown atom ',I 3 continue CLOSE(35) RETURN END c SUBROUTINE READX(R,IT,NAT,IBT) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8(A-H,O-Z) DIMENSION R(*),IT(*),IBT(*) READ(35,*) READ(35,*)NAT DO 3 I=1,NAT 3 READ(35,*)IT(I),(R(3*(I-1)+j),j=1,3),(IBT(7*(I-1)+j),j=1,7) CLOSE(35) RETURN END c SUBROUTINE WRITEX(R,IT,NAT) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8(A-H,O-Z) DIMENSION R(*),IT(*) WRITE(35,*) WRITE(35,*)NAT DO 3 I=1,NAT 3 WRITE(35,3500)IT(I),(R(3*(I-1)+j),j=1,3) 3500 FORMAT(I5,3F15.6,' 0 0 0 0 0 0 0 0.0') CLOSE(35) RETURN END c SUBROUTINE READFF(MX3,N,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(MX3,N) N1=1 1 N3=N1+4 IF(N3.GT.N)N3=N DO 130 LN=N1,N 130 READ(20,17)(FCAR(LN,J),J=N1,MIN(LN,N3)) N1=N1+5 IF(N3.LT.N)GOTO 1 17 FORMAT(4X,5D14.6) DO 31 I=1,N DO 31 J=I+1,N 31 FCAR(I,J)=FCAR(J,I) CLOSE(20) RETURN END c SUBROUTINE rCCT(NSOL,INDBIG,lcct) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION INDBIG(*) LOGICAL lcct INQUIRE(FILE='CCT.INP',exist=lcct) if(.not.lcct)return OPEN(2,FILE='CCT.INP',STATUS='OLD') READ(2,*)NAT if(NAT.NE.NSOL)THEN close(2) write(6,*)'NAT in solvent structure <> NAT/BIG in CCT.INP' stop ENDIF READ(2,*) READ(2,2001)(INDBIG(IAB),IAB=1,NAT) 2001 FORMAT(20I3) DO 5 IAB=1,NAT 5 INDBIG(IAB)=0 READ(2,2001)(INDBIG(IAB),IAB=1,NAT) c c check, that atoms are not defined twice in the overlaps: do ia=1,nat ismall=indbig(ia) if(ismall.ne.0)then do ja=ia+1,nat if(ismall.eq.indbig(ja))then write(6,*)'big ',ia,ja write(6,*)'small ',ismall write(6,*)'umbiqous overlap' stop endif enddo endif enddo c CLOSE(2) RETURN END subroutine getpot(r,ii,jj,kk,ll,fi,nat,q,AAXIS,lper) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) LOGICAL lper integer*4 q(*) real*8 r(*),fi(*) n4=4 qh=0.412d0 qo=-0.824d0 aa2=AAXIS/2.0d0 do 1 i=1,n4 if(i.eq.1)then xi=r(3*(ii-1)+1) yi=r(3*(ii-1)+2) zi=r(3*(ii-1)+3) endif if(i.eq.2)then xi=r(3*(jj-1)+1) yi=r(3*(jj-1)+2) zi=r(3*(jj-1)+3) endif if(i.eq.3)then xi=r(3*(kk-1)+1) yi=r(3*(kk-1)+2) zi=r(3*(kk-1)+3) endif if(i.eq.4)then xi=r(3*(ll-1)+1) yi=r(3*(ll-1)+2) zi=r(3*(ll-1)+3) endif fi(i)=0.0d0 do 2 j=1,nat if(q(j).eq.1.or.q(j).eq.8)then xj=r(3*(j-1)+1) yj=r(3*(j-1)+2) zj=r(3*(j-1)+3) xja=xj yja=yj zja=zj if(lper)then if(xja.gt.xi+aa2)xja=xja-AAXIS if(yja.gt.yi+aa2)yja=yja-AAXIS if(zja.gt.zi+aa2)zja=zja-AAXIS if(xja.lt.xi-aa2)xja=xja+AAXIS if(yja.lt.yi-aa2)yja=yja+AAXIS if(zja.lt.zi-aa2)zja=zja+AAXIS endif jh=0 jo=0 dij=dsqrt((xi-xja)**2+(yi-yja)**2+(zi-zja)**2) do 3 k=1,nat if((q(j).eq.1.and.q(k).eq.8).or.(q(j).eq.8.and.q(k).eq.1))then xk=r(3*(k-1)+1) yk=r(3*(k-1)+2) zk=r(3*(k-1)+3) dist=dsqrt((xk-xj)**2+(yk-yj)**2+(zk-zj)**2) if(dist.lt.1.2d0.and.dist.gt.0.8d0)then if(q(k).eq.1)jh=jh+1 if(q(k).eq.8)then kh=0 do 4 l=1,nat if(q(l).eq.1.and.l.ne.j)then xl=r(3*(l-1)+1) yl=r(3*(l-1)+2) zl=r(3*(l-1)+3) dkl=dsqrt((xk-xl)**2+(yk-yl)**2+(zk-zl)**2) if(dkl.lt.1.2d0.and.dkl.gt.0.8d0)kh=kh+1 endif 4 continue if(kh.eq.1)jo=jo+1 endif endif endif 3 continue if(q(j).eq.8.and.jh.eq.2)fi(i)=fi(i)+qo/dij if(q(j).eq.1.and.jo.eq.1)fi(i)=fi(i)+qh/dij endif 2 continue 1 continue return end c subroutine repair(F,i,j,delw,w0,MX3,r,pisvejc,e) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) REAL*8 F(MX3,MX3),e(3),r(*) c c C => O unit vector an=0.0d0 do 4 ix=1,3 e(ix)=r(3*(i-1)+ix)-r(3*(j-1)+ix) 4 an=an+e(ix)**2 do 3 ix=1,3 3 e(ix)=e(ix)/dsqrt(an) fold=0.0d0 do 1 ia=1,2 do 1 ix=1,3 if(ia.eq.1)then i1=i signi=16.0d0/28.0d0 else i1=j signi=-12.0d0/28.0d0 endif ii=3*(i1-1)+ix do 1 ja=1,2 do 1 jx=1,3 if(ja.eq.1)then j1=i signj=16.0d0/28.0d0 else j1=j signj=-12.0d0/28.0d0 endif jj=3*(j1-1)+jx 1 fold=fold+signi*e(ix)*F(ii,jj)*e(jx)*signj fac=pisvejc*delw/w0*2.0d0 fnew=fac*fold c do 2 ia=1,2 do 2 ix=1,3 if(ia.eq.1)then i1=i signi=1.0d0 else i1=j signi=-1.0d0 endif ii=3*(i1-1)+ix do 2 ja=1,2 do 2 jx=1,3 if(ja.eq.1)then j1=i signj=1.0d0 else j1=j signj=-1.0d0 endif jj=3*(j1-1)+jx 2 F(ii,jj)=F(ii,jj)+e(ix)*signi*fnew*signj*e(jx) return end subroutine repait(P,i,j,delw,w0,MX,pisvejc,e,d0,delD) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) REAL*8 P(MX,3,3),e(3) c c only diagonal elements of P uo=0.0d0 do 1 ix=1,3 uo=uo-e(ix)*P(j,ix,ix)*12.0d0/28.0d0 1 uo=uo+e(ix)*P(i,ix,ix)*16.0d0/28.0d0 un=uo*pisvejc*(delD/d0+delw/w0)/2.0d0 c do 2 ix=1,3 do 2 jj=1,3 P(i,ix,jj)=P(i,ix,jj)+e(ix)*un 2 P(j,ix,jj)=P(j,ix,jj)-e(ix)*un return end c subroutine findag(NAT,R,IT,iag,jag,kag,lag,ig,NG0) c find all mide groups in R IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) REAL*8 R(*) INTEGER*4 IT(*),iag(NG0),jag(NG0),kag(NG0),lag(NG0) c ig=0 do 1 is=1,NAT if(IT(is).eq.8)then xi=R(3*(is-1)+1) yi=R(3*(is-1)+2) zi=R(3*(is-1)+3) do 2 js=1,NAT if(IT(js).eq.6)then xj=R(3*(js-1)+1) yj=R(3*(js-1)+2) zj=R(3*(js-1)+3) d=dsqrt((xi-xj)**2+(yi-yj)**2+(zi-zj)**2) if(d.lt.1.35d0.and.d.gt.1.15d0)then do 3 ks=1,NAT if(ks.ne.js.and.IT(ks).eq.6)then xk=R(3*(ks-1)+1) yk=R(3*(ks-1)+2) zk=R(3*(ks-1)+3) d=dsqrt((xk-xj)**2+(yk-yj)**2+(zk-zj)**2) if(d.lt.1.65d0.and.d.gt.1.35d0)then do 4 ls=1,NAT if(IT(ls).eq.7)then xl=R(3*(ls-1)+1) yl=R(3*(ls-1)+2) zl=R(3*(ls-1)+3) d=dsqrt((xl-xj)**2+(yl-yj)**2+(zl-zj)**2) if(d.lt.1.48d0.and.d.gt.1.25d0)then ig=ig+1 if(ig.gt.NG0)then write(*,*)'Too many amide groups' stop endif iag(ig)=is jag(ig)=js kag(ig)=ks lag(ig)=ls endif endif 4 continue endif endif 3 continue endif endif 2 continue endif 1 continue return end subroutine ass(iags,jags,kags,lags,iagv,jagv,kagv,lagv, 1ind,il,im,RVAC,RSOL,igs,igv) c k j i l c C (C = O ) N IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) REAL*8 RVAC(*),RSOL(*) INTEGER*4 iagv(*),jagv(*),kagv(*),lagv(*),ind(*) INTEGER*4 iags(*),jags(*),kags(*),lags(*),il(*),im(*) c do 1 ig=1,igv c c find best groups to ig,igp,(igpp) il(1)=iagv(ig) il(2)=jagv(ig) il(3)=kagv(ig) il(4)=lagv(ig) natfit=4 do 11 io=1,igv if(io.ne.ig.and.natfit.le.8)then il(natfit+1)=iagv(io) il(natfit+2)=jagv(io) il(natfit+3)=kagv(io) il(natfit+4)=lagv(io) natfit=natfit+4 endif 11 continue c c try all combinations of groups in solvent del0=0.0d0 jg0=0 it=0 c c natfit= 4: ###################################### if(natfit.eq.4)then do 2 jg=1,igs im(1)=iags(jg) im(2)=jags(jg) im(3)=kags(jg) im(4)=lags(jg) call este(RVAC,RSOL,natfit,il,im,del) it=it+1 if(it.eq.1)then jg0=1 del0=del endif if(del.lt.del0)then del0=del jg0=jg endif 2 continue endif c natfit= 8: ###################################### if(natfit.eq.8)then do 298 jg=1,igs im(1)=iags(jg) im(2)=jags(jg) im(3)=kags(jg) im(4)=lags(jg) do 218 jo=1,igs if(jo.ne.jg)then im(5)=iags(jo) im(6)=jags(jo) im(7)=kags(jo) im(8)=lags(jo) it=it+1 call este(RVAC,RSOL,natfit,il,im,del) if(it.eq.1)then jg0=1 del0=del endif if(del.lt.del0)then del0=del jg0=jg endif endif 218 continue 298 continue endif c natfit=12: ###################################### if(natfit.eq.12)then do 29 jg=1,igs im(1)=iags(jg) im(2)=jags(jg) im(3)=kags(jg) im(4)=lags(jg) do 21 jo=1,igs if(jo.ne.jg)then im(5)=iags(jo) im(6)=jags(jo) im(7)=kags(jo) im(8)=lags(jo) do 22 ko=1,igs if(ko.ne.jg.and.ko.ne.jo)then im(9) =iags(ko) im(10)=jags(ko) im(11)=kags(ko) im(12)=lags(ko) it=it+1 call este(RVAC,RSOL,natfit,il,im,del) if(it.eq.1)then jg0=1 del0=del endif if(del.lt.del0)then del0=del jg0=jg endif endif 22 continue endif 21 continue 29 continue endif c ###################################### write(6,6000)ig,jg0,natfit 6000 format(i4,' assigned to ',i4,';',i4,' atom fit') ind(iags(jg0))=iagv(ig) ind(jags(jg0))=jagv(ig) ind(kags(jg0))=kagv(ig) 1 ind(lags(jg0))=lagv(ig) return end subroutine este(RVAC,RSOL,natfit,il,im,del) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (MXB=12) REAL*8 RVAC(*),RSOL(*),C1(3),C2(3) INTEGER*4 il(*),im(*) COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),RTOL,NAT,ITER DO 23 IX=1,3 C1(IX)=0.0d0 C2(IX)=0.0d0 do 23 ia=1,natfit C1(IX)=C1(IX)+RVAC(3*(il(ia)-1)+IX) 23 C2(IX)=C2(IX)+RSOL(3*(im(ia)-1)+IX) DO 24 IX=1,3 C1(IX)=C1(IX)/dble(natfit) 24 C2(IX)=C2(IX)/dble(natfit) DO 35 ISAT=1,natfit DO 35 J=1,3 XS(ISAT,J)=RVAC(3*(il(ISAT)-1)+J)-C1(J) XB(ISAT,J)=RSOL(3*(im(ISAT)-1)+J)-C2(J) 35 continue NAT=natfit CALL DOMATRIX(IERR) del=0.0d0 DO 1 I=1,NAT TX=A(1,1)*XS(I,1)+A(1,2)*XS(I,2)+A(1,3)*XS(I,3) TY=A(2,1)*XS(I,1)+A(2,2)*XS(I,2)+A(2,3)*XS(I,3) TZ=A(3,1)*XS(I,1)+A(3,2)*XS(I,2)+A(3,3)*XS(I,3) 1 del=del+(TX-XB(I,1))**2+(TY-XB(I,2))**2+(TZ-XB(I,3))**2 IF(IERR.EQ.1)then write(6,*)'Program cannot continue' stop ENDIF RETURN END c SUBROUTINE READTEN(MX,P,A,AJ,PV,NOPT,NQ) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION P(MX,3,3),A(MX,3,3),AJ(MX,3,3),PV(MX,3,3) OPEN(15,FILE='VAC.TEN') READ (15,*) NAT,NOPT,NQ DO 10 L=1,NAT DO 10 J=1,3 10 READ (15,*) (P(L,J,I),I=1,3) DO 220 L=1,NAT DO 220 J=1,3 220 READ (15,*) (A(L,J,I),I=1,3) DO 230 L=1,NAT DO 230 J=1,3 230 READ (15,*) (AJ(L,J,I),I=1,3) DO 100 L=1,NAT DO 100 J=1,3 100 READ (15,*) (PV(L,J,I),I=1,3) CLOSE(15) WRITE(*,*)' Tensors read' RETURN END C SUBROUTINE WRITETEN(MX,P,A,AJ,PV,NOPT,NQ,NAT) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION P(MX,3,3),A(MX,3,3),AJ(MX,3,3),PV(MX,3,3) OPEN(15,FILE='FILE.TEN') WRITE(15,1500) NAT,NOPT,NQ 1500 FORMAT(3I5) DO 10 L=1,NAT DO 10 J=1,3 10 WRITE(15,1501) (P(L,J,I),I=1,3),L 1501 FORMAT(3F14.8,I5) DO 220 L=1,NAT DO 220 J=1,3 220 WRITE(15,1501) (A(L,J,I),I=1,3),L DO 230 L=1,NAT DO 230 J=1,3 230 WRITE(15,1501) (AJ(L,J,I),I=1,3),L DO 100 L=1,NAT DO 100 J=1,3 100 WRITE(15,1501) (PV(L,J,I),I=1,3),L CLOSE(15) WRITE(*,*)' FILE.TEN written' RETURN END C SUBROUTINE DOMATRIX(IERR) C This is a Fortran version of Petr Malon's subroutine Amoeba IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (MXB=12) DIMENSION ANG(3),P(4,3),Y(4),PBAR(3),PR(3),PRR(3) COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),RTOL,NAT,ITER IERR=0 C STARTING VALUES FOR THE ITERATION: ANG(1)=0.2d0 ANG(2)=0.0d0 ANG(3)=0.0d0 DO 1 I=1,3 1 P(1,I)=ANG(I) Y(1)=FU(ANG) ANG(2)=ANG(1) ANG(1)=0.0d0 DO 2 I=1,3 2 P(2,I)=ANG(I) Y(2)=FU(ANG) ANG(3)=ANG(2) ANG(2)=0.0d0 DO 3 I=1,3 3 P(3,I)=ANG(I) Y(3)=FU(ANG) ANG(1)=ANG(3) ANG(2)=ANG(3) ANG(3)=0.0d0 DO 4 I=1,3 4 P(4,I)=ANG(I) Y(4)=FU(ANG) FTOL=0.0000001d0 ITMAX=1500 ITER=0 99999 ILO=1 IF(Y(1).GT.Y(2))THEN IHI=1 INHI=2 ELSE IHI=2 INHI=1 ENDIF DO 5 I=1,4 IF(Y(I).LT.Y(ILO))ILO=I IF(Y(I).GT.Y(IHI))THEN INHI=IHI IHI=I ELSE IF(Y(I).GT.Y(INHI))THEN IF(I.NE.IHI)INHI=I ENDIF ENDIF 5 CONTINUE IF(ABS(Y(IHI))+ABS(Y(ILO)).LT.0.00000000001d0)THEN RTOL=0.0d0 ELSE RTOL=ABS(Y(IHI)-Y(ILO))/(ABS(Y(IHI))+ABS(Y(ILO)))*2.0d0 ENDIF IF(RTOL.LT.FTOL)RETURN IF(ITER.EQ.ITMAX)THEN WRITE(3,*)' Rotation has not converged !' IERR=1 RETURN ENDIF ITER=ITER+1 DO 55 I=1,3 55 PBAR(I)=0.0d0 DO 6 I=1,4 IF(I.NE.IHI)THEN DO 7 J=1,3 7 PBAR(J)=PBAR(J)+P(I,J) ENDIF 6 CONTINUE DO 8 J=1,3 PBAR(J)=PBAR(J)/3.0d0 8 PR(J)=2.0d0*PBAR(J)-P(IHI,J) YPR=FU(PR) IF(YPR.LE.Y(ILO))THEN DO 9 J=1,3 9 PRR(J)=2.0d0*PR(J)-PBAR(J) YPRR=FU(PRR) IF(YPRR.LT.Y(ILO))THEN DO 10 J=1,3 10 P(IHI,J)=PRR(J) Y(IHI)=YPRR ELSE DO 11 J=1,3 11 P(IHI,J)=PR(J) Y(IHI)=YPR ENDIF ELSE IF(YPR.GE.Y(INHI))THEN IF(YPR.LT.Y(IHI))THEN DO 12 J=1,3 12 P(IHI,J)=PR(J) Y(IHI)=YPR ENDIF DO 13 J=1,3 13 PRR(J)=0.5d0*P(IHI,J)+0.5d0*PBAR(J) YPRR=FU(PRR) IF(YPRR.LT.Y(IHI))THEN DO 14 J=1,3 14 P(IHI,J)=PRR(J) Y(IHI)=YPRR ELSE DO 15 I=1,4 IF(I.NE.ILO)THEN DO 16 J=1,3 PR(J)=0.5d0*(P(I,J)+P(ILO,J)) 16 P(I,J)=PR(J) Y(I)=FU(PR) ENDIF 15 CONTINUE ENDIF ELSE DO 17 J=1,3 17 P(IHI,J)=PR(J) Y(IHI)=YPR ENDIF ENDIF GOTO 99999 END FUNCTION FU(ANG) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (MXB=12) COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),RTOL,NAT,ITER DIMENSION ANG(3) S1=SIN(ANG(1)) S2=SIN(ANG(2)) S3=SIN(ANG(3)) C1=COS(ANG(1)) C2=COS(ANG(2)) C3=COS(ANG(3)) C ang1 = theta C ang2 = phi C ang3 = cappa A(1,1)=C1*C2*C3-S2*S3 A(1,2)=C1*S2*C3+C2*S3 A(1,3)=-S1*C3 A(2,1)=-C1*C2*S3-S2*C3 A(2,2)=-C1*S2*S3+C2*C3 A(2,3)=S1*S3 A(3,1)=S1*C2 A(3,2)=S1*S2 A(3,3)=C1 R=0.0d0 DO 1 I=1,NAT TX=A(1,1)*XS(I,1)+A(1,2)*XS(I,2)+A(1,3)*XS(I,3) TY=A(2,1)*XS(I,1)+A(2,2)*XS(I,2)+A(2,3)*XS(I,3) TZ=A(3,1)*XS(I,1)+A(3,2)*XS(I,2)+A(3,3)*XS(I,3) 1 R=R+(TX-XB(I,1))*(TX-XB(I,1))+ 1 (TY-XB(I,2))*(TY-XB(I,2))+(TZ-XB(I,3))*(TZ-XB(I,3)) FU=R RETURN END C