PROGRAM hypep_ttt C empirical correction to vacuum peptide FF c on force field belonging to atoms 1...N4 c based on potential on atoms 1..N6 c nat0 - maximum number of atoms c m0 - max number of clusters c nd0 - max dimension of the inversion problem IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (nat0=800,MX3=3*nat0,MX=nat0,nd0=30) PARAMETER (nat0big=3000) DIMENSION ZIM(MX),ZM(MX3),FCAR(MX3,MX3), 2SMAT(MX3,MX3),Q(MX3),fcars(MX3,MX3),fi(nat0),r(3*nat0big), 3ity(nat0big),i4(nat0),qq(MX3),b(nat0,nd0),ind(nat0big), 4ityv(nat0),rv(MX3),qqdiff(MX3),bx(nat0,3,nd0,3),E1(MX3), 4P(nat0,3,3),aat(nat0,3,3),bp(nat0,3,nd0,3,3), 5bxm(nat0,3,nd0,3),e0(3,3),alpha(MX3,3,3),AA(MX3,3,3,3), 6bxl(nat0,3,nd0,3),G(MX3,3,3),bpk(nat0,3,nd0,3,3), 7bpm(nat0,3,nd0,3,3),bpl(nat0,3,nd0,3,3),fi_ref(NAT0) logical lex,lper,ltinker,lfdp(3) character*80 ts,name1 character*2 at DATA IERR,n6/0,6/ DATA qh,qo/0.42d0,-0.84d0/ C WRITE(6,60000) 60000 FORMAT(' PEPTIDE FORCE FIELD SOLVENT PERTURBATION ',/,/, 1 ' Input: BIJ.SCR (fitting coeffs.) ',/, 1 ' SMALL.FC (vacuum ff) ',/, 1 ' SMALL.X (vacuum geo) ',/, 1 ' CCT.INP (atom map) ',/, 1 ' BIG.X ',/, 1 ' HYPEP.OPT (optional) ',/, 1 ' Output: FILE.FC ',/, 1 ' FILE.TEN (vac, still the same) ',/, 1 ' FILE.X ',/,/) C open(55,file='BIJ.SCR',form='unformatted',status='old') read(55)n4,n6,(zim(i),i=1,n4),(lfdp(i),i=1,3) if(lfdp(1))read(55)((b(i,j),i=1,3*n4),j=1,n6) if(lfdp(2))read(55) 1((((bx(i,i3,j,ix),i=1,n4),i3=1,3),j=1,n6),ix=1,3) if(lfdp(3))read(55) 1(((((bp(i,i3,j,ix,iz),i=1,n4),i3=1,3),j=1,n6),ix=1,3),iz=1,3) read(55)n6,(fi_ref(i),i=1,n6) close(55) DO 220 i=1,n4 DO 220 j=1,3 220 ZM(3*(i-1)+j)=1.0d0/DSQRT(zim(i)) write(6,*)'Fitting constants read from BIJ.SCR.' write(6,*)lfdp c OPEN(20,FILE='SMALL.X',STATUS='OLD') read(20,*) read(20,*)natv if(natv.gt.nat0)call report('Too many atoms') do 41 ia=1,natv 41 read(20,*)ityv(ia),(rv(3*(ia-1)+ix),ix=1,3) close(20) c OPEN(20,FILE='FILE.X') write(20,*)'Vacuum peptide geometry' write(20,*)natv do 411 ia=1,natv 411 write(20,200)ityv(ia),(rv(3*(ia-1)+ix),ix=1,3) 200 format(i4,3f15.8,' 0 0 0 0 0 0 0 0.0') close(20) write(6,*)'Vacuum geometry rewritten into FILE.X' lper=.false. ltinker=.false. name1='BIG.X' inquire(file='HYPEP.OPT',exist=lex) if(lex)then open(55,file='HYPEP.OPT') 1001 read(55,2000,end=1002,err=1002)ts 2000 format(a80) if(ts(1:4).eq.'NAME')read(55,2000)name1 if(ts(1:4).eq.'TINK')then ltinker=.true. write(6,*)'Tinker format for BIG' name1='tinker.xyz' endif if(ts(1:4).eq.'PERI')then lper=.true. read(55,*)DPCX,DPCY,DPCZ dx2=DPCX/2.0d0 dy2=DPCY/2.0d0 dz2=DPCZ/2.0d0 write(6,*)' Periodic boundary conditions' endif c identify proline (changed) atoms in the MD file if(ts(1:4).eq.'CHAN')then read(55,*)n6 if(n6.gt.nat0)call report('too many changed atoms') read(55,*)(i4(i),i=1,n6) endif goto 1001 1002 close(55) write(6,*)' HYPEP.OPT found' endif c OPEN(20,FILE=name1,STATUS='OLD') if(.not.ltinker)read(20,*) read(20,*)nat if(nat.gt.nat0big)call report('Too many atoms') if(ltinker)then do 4141 ia=1,nat read(20,2020)at,(r(3*(ia-1)+ix),ix=1,3) 2020 format(6x,2x,a2,1x,3f12.6) ity(ia)=0 if(at(1:1).eq.'H')ity(ia)=1 if(at(1:1).eq.'C')ity(ia)=6 if(at(1:1).eq.'N')ity(ia)=7 if(at(1:1).eq.'O')ity(ia)=8 4141 if(ity(ia).eq.0)write(6,*)'Unknown atom ',ia,at else do 414 ia=1,nat 414 read(20,*)ity(ia),(r(3*(ia-1)+ix),ix=1,3) endif close(20) write(6,*)'Geometry read for BIG.' c OPEN(20,FILE='CCT.INP',STATUS='OLD') read(20,*) read(20,*) read(20,201)(ind(i),i=1,nat) read(20,201)(ind(i),i=1,nat) 201 format(20i3) close(20) write(6,*)' Atom map read from CCT.INP, nat= ',nat c C Read-in vacuum Cartesian FF if(lfdp(1))then OPEN(20,FILE='SMALL.FC',STATUS='OLD') CALL READFF(MX3,3*natv,FCAR) CLOSE(20) endif c C Read-in vacuum Cartesian dipole derivatives if(lfdp(2))then OPEN(15,FILE='SMALL.TEN',STATUS='OLD') CALL READTEN(nat0,P,aat,natv) CLOSE(15) endif c C Read-in vacuum Cartesian polarizability derivatives if(lfdp(3))then OPEN(2,FILE='SMALL.TTT',STATUS='OLD') call READTTT(MX3,natv,alpha,G,AA) CLOSE(2) endif c c identify proline (changed) atoms c c potential on "changed" atoms: do 4 i=1,n6 4 fi(i)=0.0d0 c c find waters and their potentials do 6 jj=1,nat if(ity(jj).eq.8)then j1=3*(jj-1) xo=r(j1+1) yo=r(j1+2) zo=r(j1+3) nh=0 do 7 kk=1,nat if(ity(kk).eq.1)then k1=3*(kk-1) xh=r(k1+1) yh=r(k1+2) zh=r(k1+3) d2=(xo-xh)**2+(yo-yh)**2+(zo-zh)**2 if(d2.gt.0.64d0.and.d2.lt.1.44d0)then nh=nh+1 if(nh.gt.2)call report('Too many Hs around one O') if(nh.eq.1)then x1=xh y1=yh z1=zh else x2=xh y2=yh z2=zh endif endif endif 7 continue if(nh.eq.2)then do ii=1,n6 call add(i4(ii) ,ii,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo, 1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ) enddo endif endif 6 continue write(6,*)' Potentials:' write(6,1603)(fi(ii),ii=1,N6) 1603 format(8f9.3) do ii=1,N6 fi(ii)=fi(ii)-fi_ref(ii) enddo write(6,*)' Reference:' write(6,1603)(fi_ref(ii),ii=1,N6) write(6,*)' Difference from the reference:' write(6,1603)(fi(ii),ii=1,N6) write(6,*) write(6,*)'group in MD:' write(6,67900)(i4(ii),ii=1,n4) 67900 format(20i4) write(6,*)'group in SMALL molecule:' write(6,67900)(ind(i4(ii)),ii=1,n4) c if(lfdp(1))then c c rewrite vac FF into smaller matrix do 31 ii=1,n4 ir=ind(i4(ii)) if(ir.eq.0)call report('Undefined overlap') do 31 ix=1,3 iin=3*(ii-1)+ix i1 =3*(ir-1)+ix do 31 jj=1,n4 ja=ind(i4(jj)) do 31 jx=1,3 jin=3*(jj-1)+jx j1 =3*(ja-1)+jx 31 fcars(iin,jin)=FCAR(i1,j1) c calculate eigenvalues IMOD=2 NA=3*N4 DO 214 i1=1,NA DO 214 j1=1,NA 214 fcars(i1,j1)=fcars(i1,j1)*zm(j1)*zm(i1) CALL TRED12(MX3,NA,fcars,SMAT,Q,IMOD,IERR,E1) IF(IERR.NE.0)call report(' Diagonalization error') DO 310 Iz=1,NA t=Q(Iz) qq(Iz)=t Q(Iz)=SQRT(ABS(t))*1302.828d0 310 if(t.lt.0.0d0)Q(Iz)=-Q(Iz) write(6,*)' Eigenvalues:' write(6,1604)(Q(II),II=1,NA) 1604 format(8f9.1) c c repaired eigenvalues - only difference to vacuum! do 100 Iz=1,NA qqdiff(Iz)=0.0d0 do 100 i=1,n6 100 qqdiff(Iz)=qqdiff(Iz)+b(Iz,i)*fi(i) c DO 312 Iz=1,NA t=qq(Iz)+qqdiff(Iz) Q(Iz)=SQRT(ABS(t))*1302.828d0 312 if(t.lt.0.0d0)Q(Iz)=-Q(Iz) write(6,*)' New eigenvalues:' write(6,1604)(Q(II),II=1,NA) c c repaired force constants do 101 i1=1,NA do 101 j1=1,NA sum=0.0d0 do 102 Iz=1,NA 102 sum=sum+SMAT(i1,Iz)*qqdiff(Iz)*SMAT(j1,Iz) 101 fcars(i1,j1)=fcars(i1,j1)+sum c DO 215 i1=1,NA DO 215 j1=1,NA 215 fcars(i1,j1)=fcars(i1,j1)/(zm(j1)*zm(i1)) c c transcript back do 311 ii=1,n4 ir=ind(i4(ii)) do 311 ix=1,3 iin=3*(ii-1)+ix i1 =3*(ir-1)+ix do 311 jj=1,n4 ja=ind(i4(jj)) do 311 jx=1,3 jin=3*(jj-1)+jx j1 =3*(ja-1)+jx 311 FCAR(i1,j1)=fcars(iin,jin) endif c c local coordinate system do 50 ix=1,3 e0(1,ix)=rv(3*(ind(i4(4))-1)+ix)-rv(3*(ind(i4(3))-1)+ix) 50 e0(3,ix)=rv(3*(ind(i4(2))-1)+ix)-rv(3*(ind(i4(3))-1)+ix) call norm(e0,1) e0(2,1)=e0(1,2)*e0(3,3)-e0(1,3)*e0(3,2) e0(2,2)=e0(1,3)*e0(3,1)-e0(1,1)*e0(3,3) e0(2,3)=e0(1,1)*e0(3,2)-e0(1,2)*e0(3,1) call norm(e0,2) e0(3,1)=e0(1,2)*e0(2,3)-e0(1,3)*e0(2,2) e0(3,2)=e0(1,3)*e0(2,1)-e0(1,1)*e0(2,3) e0(3,3)=e0(1,1)*e0(2,2)-e0(1,2)*e0(2,1) c if(lfdp(2))then c c transform coefficients into laboratory system do 551 j=1,n6 do 551 i=1,n4 c c first index do 5511 ix=1,3 do 5511 iz=1,3 sum=0.0d0 do 5512 l=1,3 5512 sum=sum+bx(i,l,j,ix)*e0(l,iz) 5511 bxm(i,iz,j,ix)=sum c c second index do 551 ix=1,3 do 551 iz=1,3 sum=0.0d0 do 5513 l=1,3 5513 sum=sum+bxm(i,ix,j,l)*e0(l,iz) 551 bxl(i,ix,j,iz)=sum c c repair APT do 1022 i=1,n4 iav=ind(i4(i)) do 1022 ix=1,3 do 1022 iy=1,3 do 1022 j=1,n6 1022 P(iav,iy,ix)=P(iav,iy,ix)+bxl(i,iy,j,ix)*fi(j) endif if(lfdp(3))then c c transform coefficients into laboratory system do 651 j=1,n6 do 651 i=1,n4 c c first index do 6511 ix=1,3 do 6511 iz=1,3 do 6511 iy=1,3 sum=0.0d0 do 6512 l=1,3 6512 sum=sum+bp(i,l,j,ix,iy)*e0(l,iz) 6511 bpm(i,iz,j,ix,iy)=sum c c second index do 652 ix=1,3 do 652 iz=1,3 do 652 iy=1,3 sum=0.0d0 do 6513 l=1,3 6513 sum=sum+bpm(i,ix,j,l,iy)*e0(l,iz) 652 bpl(i,ix,j,iz,iy)=sum c c third index do 651 ix=1,3 do 651 iz=1,3 do 651 iy=1,3 sum=0.0d0 do 6514 l=1,3 6514 sum=sum+bpl(i,ix,j,iz,l)*e0(l,iz) 651 bpk(i,ix,j,iz,iy)=sum c c repair alpha do 2022 i=1,n4 iav=ind(i4(i)) do 2022 ix=1,3 ii=ix+3*(iav-1) do 2022 iy=1,3 do 2022 iz=1,3 do 2022 j=1,n6 2022 alpha(ii,iy,iz)=alpha(ii,iy,iz)+bpk(i,ix,j,iy,iz)*fi(j) endif c c c save repaired tensors: if(lfdp(1))call WRITEFF(MX3,3*natv,FCAR) if(lfdp(2))call WRITETEN(nat0,P,aat,natv) if(lfdp(3))call WRITETTT(MX3,natv,alpha,AA,G,alpha) c stop c end SUBROUTINE WRITETTT(N30,NAT,ALPHA,A,G,alphav) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION ALPHA(N30,3,3),G(N30,3,3),A(N30,3,3,3), 1alphav(N30,3,3) c OPEN(22,FILE='FILE.TTT') WRITE(22,2000)NAT 2000 FORMAT(' ROA tensors, cartesian derivatives',/,I4,' atoms',/, 1' The electric-dipolar electric-dipolar polarizability:',/, 2' Atom/x jx jy jz') DO 1 I=1,3 WRITE(22,2002)I 2002 FORMAT(' Alpha(',I1,',J):') DO 1 L=1,NAT DO 1 IX=1,3 II=3*(L-1)+IX 1 WRITE(22,2001)L,IX,(ALPHA(II,I,J),J=1,3) 2001 FORMAT(I5,1H ,I1,3F13.7) WRITE(22,2004) 2004 FORMAT(' The electric dipole magnetic dipole polarizability:',/, 1 ' Atom/x jx(Bx) jy(By) jz(Bz)') DO 2 I=1,3 WRITE(22,2003)I 2003 FORMAT(' G(',I1,',J):') DO 2 L=1,NAT DO 2 IX=1,3 II=3*(L-1)+IX c last inner index magnetic 2 WRITE(22,2001)L,IX,(G(II,I,J),J=1,3) WRITE(22,2005) 2005 FORMAT(' The electric dipole electric quadrupole polarizability:', 2/, ' Atom/x kx ky kz') DO 3 I=1,3 DO 3 J=1,3 WRITE(22,2006)I,J 2006 FORMAT(' A(',I1,'(u),',I1,',K):') DO 3 L=1,NAT DO 3 IX=1,3 II=3*(L-1)+IX 3 WRITE(22,2001)L,IX,(A(II,I,J,K),K=1,3) WRITE(22,2009) 2009 FORMAT('The velocity form of alpha',/, 1' Atom/x vx vy vz') DO 11 I=1,3 WRITE(22,2012)I 2012 FORMAT(' AlphaV(',I1,',J):') DO 11 L=1,NAT DO 11 IX=1,3 II=3*(L-1)+IX 11 WRITE(22,2001)L,IX,(ALPHAV(II,I,J),J=1,3) CLOSE(22) write(6,*)'File FILE.TTT written' RETURN END subroutine norm(v,i) real*8 v(3,3),an integer*4 i an=dsqrt(v(i,1)**2+v(i,2)**2+v(i,3)**2) v(i,1)=v(i,1)/an v(i,2)=v(i,2)/an v(i,3)=v(i,3)/an return end subroutine report(s) character*(*) s write(6,*)s stop end subroutine add(ii,ic,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo, 1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ) real*8 r(*),xo,yo,zo,x1,y1,z1,x2,y2,z2,fi(*),x,y,z,ro,r1,r2,qh,qo, 1dx2,dy2,dz2,DPCX,DPCY,DPCZ integer*4 ii,ic logical lper c x=r(3*(ii-1)+1) y=r(3*(ii-1)+2) z=r(3*(ii-1)+3) if(lper)then if(x-xo.gt.dx2.and.x-x1.gt.dx2.and.x-x2.gt.dx2)x=x-DPCX if(y-yo.gt.dy2.and.y-y1.gt.dy2.and.y-y2.gt.dy2)y=y-DPCY if(z-zo.gt.dz2.and.z-z1.gt.dz2.and.z-z2.gt.dz2)z=z-DPCZ if(xo-x.gt.dx2.and.x1-x.gt.dx2.and.x2-x.gt.dx2)x=x+DPCX if(yo-y.gt.dy2.and.y1-y.gt.dy2.and.y2-y.gt.dy2)y=y+DPCY if(zo-z.gt.dz2.and.z1-z.gt.dz2.and.z2-z.gt.dz2)z=z+DPCZ endif ro=dsqrt((xo-x)**2+(yo-y)**2+(zo-z)**2) r1=dsqrt((x1-x)**2+(y1-y)**2+(z1-z)**2) r2=dsqrt((x2-x)**2+(y2-y)**2+(z2-z)**2) c fi(ic)=fi(ic)+qh/r1+qh/r2+qo/ro return end C SUBROUTINE WRITETEN(N0,P,A,NAT) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION P(N0,3,3),A(N0,3,3) OPEN(15,FILE='FILE.TEN') WRITE(15,1500) NAT,NAT-6,0 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 Z0=0.0d0 DO 230 L=1,NAT DO 230 J=1,3 230 WRITE(15,1501) (Z0,I=1,3),L DO 100 L=1,NAT DO 100 J=1,3 100 WRITE(15,1501) (P(L,J,I),I=1,3),L CLOSE(15) write(6,*)'File FILE.TEN written' RETURN END c SUBROUTINE TRED12(MX3,N,A,Z,D,IEIG,IERR,E) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION A(MX3,N),Z(MX3,N),D(MX3),E(MX3) IF (IEIG .LT. 0) GO TO 110 DO 100 I = 1, N C DO 100 J = 1, I Z(I,J) = A(I,J) 100 CONTINUE C 110 ICODE = IABS(IEIG) IF (N .EQ. 1) GO TO 320 C :::::::::: FOR I=N STEP -1 UNTIL 2 DO -- :::::::::: DO 300 II = 2, N I = N + 2 - II L = I - 1 H = 0.0D0 SCALE = 0.0D0 IF (L .LT. 2) GO TO 130 C :::::::::: SCALE ROW (ALGOL TOL THEN NOT NEEDED) :::::::::: DO 120 K = 1, L 120 SCALE = SCALE + DABS(Z(I,K)) C IF (DABS(SCALE) .GT. 1.0D-10) GO TO 140 130 E(I) = Z(I,L) GO TO 290 C 140 DO 150 K = 1, L F = Z(I,K) / SCALE Z(I,K) = F H = H + F * F 150 CONTINUE C G = -DSIGN(DSQRT(H),F) E(I) = SCALE * G H = H - F * G Z(I,L) = F - G F = 0.0D0 C DO 240 J = 1, L IF (ICODE .EQ. 2) Z(J,I) = Z(I,J) / (SCALE * H) G = 0.0D0 C :::::::::: FORM ELEMENT OF A*U :::::::::: DO 180 K = 1, J 180 G = G + Z(J,K) * Z(I,K) C JP1 = J + 1 IF (L .LT. JP1) GO TO 220 C DO 200 K = JP1, L 200 G = G + Z(K,J) * Z(I,K) C :::::::::: FORM ELEMENT OF P :::::::::: 220 E(J) = G / H F = F + E(J) * Z(I,J) 240 CONTINUE C HH = F / (H + H) C :::::::::: FORM REDUCED A :::::::::: DO 260 J = 1, L F = Z(I,J) G = E(J) - HH * F E(J) = G C DO 260 K = 1, J Z(J,K) = Z(J,K) - F * E(K) - G * Z(I,K) 260 CONTINUE C DO 280 K = 1, L 280 Z(I,K) = SCALE * Z(I,K) C 290 D(I) = H 300 CONTINUE C 320 D(1) = 0.0D0 E(1) = 0.0D0 C :::::::::: ACCUMULATION OF TRANSFORMATION MATRICES :::::::::: IF (ICODE .NE. 2) GO TO 600 DO 500 I = 1, N L = I - 1 IF (DABS(D(I)) .LT. 1.0D-10) GO TO 380 C DO 360 J = 1, L G = 0.0D0 C DO 340 K = 1, L 340 G = G + Z(I,K) * Z(K,J) C DO 360 K = 1, L Z(K,J) = Z(K,J) - G * Z(K,I) 360 CONTINUE C 380 D(I) = Z(I,I) Z(I,I) = 1.0D0 IF (L .LT. 1) GO TO 500 C DO 400 J = 1, L Z(I,J) = 0.0D0 Z(J,I) = 0.0D0 400 CONTINUE C 500 CONTINUE C 620 CALL TQL12 (MX3,N,Z,D,IERR,ICODE,E) RETURN C :::::::::: ALTERNATE FINAL LOOP FOR EIGENVALUES ONLY :::::::::: 600 DO 610 I=1,N 610 D(I) = Z(I,I) GO TO 620 END C SUBROUTINE TQL12(MX3,N,Z,D,IERR,ICODE,E) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION Z(MX3,N),D(MX3),E(MX3) REAL*8 MACHEP EPS = 1.0D0 10 EPS = 0.50D0*EPS TOL1 = EPS + 1.0D0 IF((TOL1.GT.1.0D0).AND.(TOL1-EPS.EQ.1.0D0)) GO TO 10 IF(TOL1-EPS.EQ.1.0D0) EPS = EPS + EPS MACHEP = EPS C C MACHEP=16.0D0**(-13) C IERR = 0 IF (N .EQ. 1) GO TO 1001 C DO 100 I = 2, N 100 E(I-1) = E(I) C F = 0.0D0 B = 0.0D0 E(N) = 0.0D0 C DO 240 L = 1, N J = 0 H = MACHEP * (DABS(D(L)) + DABS(E(L))) IF (B .LT. H) B = H C :::::::::: LOOK FOR SMALL SUB-DIAGONAL ELEMENT :::::::::: DO 110 M = L, N IF (DABS(E(M)) .LE. B) GO TO 120 C :::::::::: E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT C THROUGH THE BOTTOM OF THE LOOP :::::::::: 110 CONTINUE C 120 IF (M .EQ. L) GO TO 220 130 IF (J .EQ. 30) GO TO 1000 J = J + 1 C :::::::::: FORM SHIFT :::::::::: L1 = L + 1 G = D(L) P = (D(L1) - G) / (2.0D0 * E(L)) R = DSQRT(P*P+1.0D0) D(L) = E(L) / (P + DSIGN(R,P)) H = G - D(L) C DO 140 I = L1, N 140 D(I) = D(I) - H C F = F + H C :::::::::: QL TRANSFORMATION :::::::::: P = D(M) C = 1.0D0 S = 0.0D0 MML = M - L C :::::::::: FOR I=M-1 STEP -1 UNTIL L DO -- :::::::::: DO 200 II = 1, MML I = M - II G = C * E(I) H = C * P IF (DABS(P) .LT. DABS(E(I))) GO TO 150 C = E(I) / P R = DSQRT(C*C+1.0D0) E(I+1) = S * P * R S = C / R C = 1.0D0 / R GO TO 160 150 C = P / E(I) R = DSQRT(C*C+1.0D0) E(I+1) = S * E(I) * R S = 1.0D0 / R C = C * S 160 P = C * D(I) - S * G D(I+1) = H + S * (C * G + S * D(I)) C :::::::::: FORM VECTOR :::::::::: IF (ICODE .NE. 2) GO TO 200 DO 180 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 180 CONTINUE C 200 CONTINUE C E(L) = S * P D(L) = C * P IF (DABS(E(L)) .GT. B) GO TO 130 220 D(L) = D(L) + F 240 CONTINUE C :::::::::: ORDER EIGENVALUES AND EIGENVECTORS :::::::::: DO 300 II = 2, N I = II - 1 K = I P = D(I) C DO 260 J = II, N IF (D(J) .LE. P) GO TO 260 K = J P = D(J) 260 CONTINUE C IF (K .EQ. I) GO TO 300 D(K) = D(I) D(I) = P C IF (ICODE .NE. 2) GO TO 300 DO 280 J = 1, N P = Z(J,I) Z(J,I) = Z(J,K) Z(J,K) = P 280 CONTINUE C 300 CONTINUE C GO TO 1001 C :::::::::: SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30 ITERATIONS :::::::::: 1000 IERR = L 1001 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) C CONST=4.359828d0/0.5291772d0/0.5291772d0 DO 3 I=1,N DO 3 J=I,N 3 FCAR(J,I)=FCAR(J,I)*CONST DO 31 I=1,N DO 31 J=I+1,N 31 FCAR(I,J)=FCAR(J,I) WRITE(*,*)' Cartesian FF read in ... ' RETURN END 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') CONST=4.359828d0/0.5291772d0**2 DO 6 I=1,N DO 6 J=1,N 6 FCAR(I,J)=FCAR(I,J)/CONST 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) write(6,*)'Repaired FF written into FILE.FC.' RETURN END SUBROUTINE READTTT(N30,NAT,ALPHA,G,A) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION ALPHA(N30,3,3),G(N30,3,3),A(N30,3,3,3) READ(2,*) READ(2,*)NAT READ(2,*) READ(2,*) DO 1 I=1,3 READ(2,*) DO 1 L=1,NAT DO 1 IX=1,3 IND=3*(L-1)+IX 1 READ(2,*)LDUM,IXDUM,(ALPHA(IND,I,J),J=1,3) READ(2,*) READ(2,*) DO 2 I=1,3 READ(2,*) DO 2 L=1,NAT DO 2 IX=1,3 IND=3*(L-1)+IX c last index magnetic, inner 2 READ(2,*)LDUM,IXDUM,(G(IND,I,J),J=1,3) READ(2,*) READ(2,*) DO 3 I=1,3 DO 3 J=1,3 READ(2,*) DO 3 L=1,NAT DO 3 IX=1,3 IND=3*(L-1)+IX 3 READ(2,*)LDUM,IXDUM,(A(IND,I,J,K),K=1,3) RETURN END SUBROUTINE READTEN(N0,P,A,NAT) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION P(N0,3,3),A(N0,3,3) READ(15,*) 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,*) DO 100 L=1,NAT DO 100 J=1,3 100 READ (15,*) RETURN END