PROGRAM hypep 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=10) 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) logical lex,lper,ltinker,lfdp(3),lter 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) close(55) open(55,file='BIJ.TXT') write(55,555)n4,n6 write(55,552)(zim(i),i=1,n4) write(55,553)(lfdp(i),i=1,3) 555 format(2i6) 552 format(8f10.1) 553 format(3l2) if(lfdp(1))then do i=1,3*n4 write(55,554)(b(i,j),j=1,n6) 554 format(8G10.2) enddo endif 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. c c if also tertiary amide group: lter=.true. 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.'TERT')read(55,*)lter 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 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' 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 amide groups in peptide iag=0 do 1 ia=1,nat if(ity(ia).eq.6)then xi=r(3*(ia-1)+1) yi=r(3*(ia-1)+2) zi=r(3*(ia-1)+3) iin=0 iih=0 iian1=0 iian2=0 iiac=0 iio=0 in=0 ih=0 ian=0 iac=0 io=0 it=0 do 2 ja=1,nat if(ja.ne.ia)then xj=r(3*(ja-1)+1) yj=r(3*(ja-1)+2) zj=r(3*(ja-1)+3) d2=(xi-xj)**2+(yi-yj)**2+(zi-zj)**2 if(d2.lt.2.7d0)then it=it+1 if(ity(ja).eq.8)then io=io+1 iio=ja endif if(ity(ja).eq.6)then iac=iac+1 iiac=ja endif if(ity(ja).eq.7)then in=in+1 iin=ja endif endif endif 2 continue c c 3 bonds: to O C N if(it.eq.3.and.io.eq.1.and.iac.eq.1.and.in.eq.1)then xoo=r(3*(iio-1)+1) yoo=r(3*(iio-1)+2) zoo=r(3*(iio-1)+3) xn=r(3*(iin-1)+1) yn=r(3*(iin-1)+2) zn=r(3*(iin-1)+3) do 3 ja=1,nat if(ja.ne.iin)then xj=r(3*(ja-1)+1) yj=r(3*(ja-1)+2) zj=r(3*(ja-1)+3) d2=(xn-xj)**2+(yn-yj)**2+(zn-zj)**2 if(d2.lt.2.7d0)then if(ity(ja).eq.1)then ih=ih+1 iih=ja endif if(ity(ja).eq.6.and.ja.ne.ia)then ian=ian+1 if(iian1.eq.0)then iian1=ja else iian2=ja endif endif endif endif 3 continue if(ih.eq.1.and.ian.eq.1)write(6,7001) 7001 format('Secondary amide:') if(ih.eq.0.and.ian.eq.2)then if(lter)then write(6,7002) 7002 format('Tertiary amide:') c determina, which carbon is at the hydrogen position:' xtt=r(3*(iian1-1)+1) ytt=r(3*(iian1-1)+2) ztt=r(3*(iian1-1)+3) d1=(xoo-xtt)**2+(yoo-ytt)**2+(zoo-ztt)**2 xtt=r(3*(iian2-1)+1) ytt=r(3*(iian2-1)+2) ztt=r(3*(iian2-1)+3) d2=(xoo-xtt)**2+(yoo-ytt)**2+(zoo-ztt)**2 if(d2.lt.d1)then item=iian1 iian1=iian2 iian2=item endif else write(6,7003) 7003 format('Tertiary amide encountered, will be skipped') endif endif if((ih.eq.1.and.ian.eq.1). 1 or.(ih.eq.0.and.ian.eq.2.and.lter))then iag=iag+1 iic=ia imap=0 if(ian.eq.1)then write(6,600)iag,iih,iin,iic,iio,iian1,iiac write(6,601)ind(iih),ind(iin),ind(iic),ind(iio),ind(iian1), 1 ind(iiac) if(ind(iih).eq.0.or.ind(iin).eq.0.or.ind(iic).eq.0.or. 1 ind(iio).eq.0.or.ind(iian1).eq.0.or.ind(iiac).eq.0)then write(6,65) 65 format('Amide group ill-mapped') imap=1 endif else write(6,600)iag,iian2,iin,iic,iio,iian1,iiac write(6,601)ind(iian2),ind(iin),ind(iic),ind(iio),ind(iian1), 1 ind(iiac) if(ind(iian2).eq.0.or.ind(iin).eq.0.or.ind(iic).eq.0.or. 1 ind(iio).eq.0.or.ind(iian1).eq.0.or.ind(iiac).eq.0)then write(6,65) imap=1 endif endif 600 format(' Amide group number ',i3,': ',6i3) 601 format(' Vacuum system ',3x,': ',6i3) if(imap.eq.0)then if(ian.eq.1)then i4(1)=iih else i4(1)=iian2 endif i4(2)=iin i4(3)=iic i4(4)=iio if(n4.eq.6)then i4(5)=iian1 i4(6)=iiac endif c c find waters and their potentials do 4 i=1,n6 4 fi(i)=0.0d0 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 if(ian.eq.1)then call add(iih ,1,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo, 1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ) else call add(iian2,1,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo, 1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ) endif call add(iin ,2,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo, 1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ) call add(iic ,3,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo, 1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ) call add(iio ,4,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo, 1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ) call add(iian1,5,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo, 1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ) call add(iiac,6,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo, 1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ) endif endif 6 continue write(6,*)' Potentials:' write(6,1603)(fi(ii),ii=1,N6) 1603 format(8f9.3) 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)then write(6,*)'n4 ii ir ',n4,ii,ir,i4(ii) call report('Undefined overlap') endif 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 c c skip "H" in tertiary amides: if(ian.eq.1.or.(i1.gt.3.and.j1.gt.3))then sum=0.0d0 do 102 Iz=1,NA 102 sum=sum+SMAT(i1,Iz)*qqdiff(Iz)*SMAT(j1,Iz) fcars(i1,j1)=fcars(i1,j1)+sum endif 101 continue 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 10221 i=1,n4 c skip "H" in tertiary amides: if(ian.eq.1.or.i.gt.1)then 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 10221 continue 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 20221 i=1,n4 c skip "H" in tertiary amides: if(ian.eq.1.or.i.gt.1)then 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 20221 continue endif c endif endif endif endif 1 continue write(6,*)iag,' amide groups found' 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