PROGRAM FITNMA C finds fitting coefficients for solvation effect c on force field belonging to atoms 1...N4 c based on potential on atoms 1..N6 IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION e0(3,3) CHARACTER*80 sname CHARACTER*80,allocatable::ffield(:),xnamel(:),tnamel(:),pnamel(:) real*8, allocatable::r(:,:),fi(:,:),fi_ref(:), 1FCAR(:,:),fcars(:,:),P(:,:,:),AAT(:,:,:), 1qqx(:,:,:,:),alpha(:,:,:),G(:,:,:),AT(:,:,:,:),qt(:,:,:), 6qpx(:,:,:,:,:),a(:,:),work(:),am(:,:),c(:),b(:,:,:),d(:,:), 1bx(:,:,:,:),bp(:,:,:,:,:),bg(:,:,:,:,:),ba(:,:,:,:,:,:), 6qgx(:,:,:,:,:),qg(:,:,:),qa(:,:,:,:),ax(:,:,:,:,:,:),ff(:,:,:) integer*4, allocatable::i4(:,:),i6(:,:),ity(:),ipiv(:) LOGICAL lfdp(3),detect,lcompl DATA qh,qo/0.42d0,-0.84d0/ C WRITE(6,60000) 60000 FORMAT(' FORCE FIELD SOLVENT PERTURBATION',/,/, 1 ' Input : FITNMA.INP ',/, 1 ' Output: BIJ.SCR ',/) C OPEN(16,FILE='FITNMA.INP',status='old') read(16,1600)sname write(6,1600)sname 1600 format(a80) read(16,*)m write(6,1601)m 1601 format(i5,' md configurations') read(16,*)n4,n6 write(6,1602)n4,n6 1602 format(i5,' perturbed atoms, ',i2,' potential sites') c c read masses, just for units of the FCs: allocate(fi_ref(n6)) c c read whether to fit fc, P, A: read(16,*)(lfdp(i),i=1,3) write(6,600)(lfdp(i),i=1,3) 600 format(' Fit force field, dipole and polarizability dervs.:',3l2) allocate(ffield(m),xnamel(m),tnamel(m),pnamel(m),i6(m,n6), 1i4(m,n4)) im=0 nmax=0 do 1 i=1,m c FC, dipole derivaties, polarizability derivatives, coords: if(lfdp(1))read(16,1600)ffield(i) if(lfdp(2))read(16,1600)tnamel(i) if(lfdp(3))read(16,1600)pnamel(i) read(16,1600)xnamel(i) call rgn(xnamel(i),nat) read(16,*)(i4(i,j),j=1,n4),(i6(i,j),j=1,n6) if(nat.gt.nmax)then nmax=nat im=i endif 1 continue close(16) write(6,*)'Maximal number of atoms',nmax,', cluster',im c ========================================================= allocate(ity(nmax),fi(m,nmax), 1fcars(3*n4,3*n4),P(nmax,3,3), 1AAT(nmax,3,3),qqx(m,nmax,3,3),ff(m,3*n4,3*n4), 6qpx(m,2*n4,3,3,3),qgx(m,2*n4,3,3,3),ax(m,2*n4,3,3,3,3)) fi=0.0d0 do 2 i=1,m call rgn(xnamel(i),nat) allocate(alpha(6*nat,3,3),AT(6*nat,3,3,3),G(6*nat,3,3), 1 qt(6*nat,3,3),qg(6*nat,3,3),qa(6*nat,3,3,3),r(nat,3), 1 FCAR(3*nat,3*nat)) c read geometry call rg(xnamel(i),nat,ity,r,1) c local coordinate system, first three atoms i4: call lcs(e0,r,nmax,i4,n4,m,i) c c identify waters: do 6 jj=1,nat if(ity(jj).eq.8)then xo=r(jj,1) yo=r(jj,2) zo=r(jj,3) nh=0 do 7 kk=1,nat if(ity(kk).eq.1)then xh=r(kk,1) yh=r(kk,2) zh=r(kk,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 c add to potentials on atoms 1..N6 do 5 ii=1,N6 x=r(i6(i,ii),1) y=r(i6(i,ii),2) z=r(i6(i,ii),3) 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) 5 fi(i,ii)=fi(i,ii)+qh/r1+qh/r2+qo/ro endif endif 6 continue write(6,*)' Potentials:' write(6,1607)(fi(i,ii),ii=1,N6) 1607 format(8f10.3) c c first cluster as a reference: if(i.eq.1)then do 8 ii=1,N6 8 fi_ref(ii)=fi(i,ii) else write(6,*)' Difference vs reference (first cluster):' do 9 ii=1,N6 9 fi(i,ii)=fi(i,ii)-fi_ref(ii) write(6,1607)(fi(i,ii),ii=1,N6) endif c C Read-in cartesian FF if(lfdp(1))then CALL READFF(3*nat,FCAR,ffield(i)) c rewrite into small matrix do 3 ii=1,n4 ia=i4(i,ii) do 3 ix=1,3 iin=3*(ii-1)+ix i1 =3*(ia-1)+ix do 3 jj=1,n4 ja=i4(i,jj) do 3 jx=1,3 jin=3*(jj-1)+jx j1 =3*(ja-1)+jx 3 fcars(iin,jin)=FCAR(i1,j1) call trf(n4,fcars,e0,ia) do 31 ii=1,n4 do 31 jj=1,n4 31 ff(i,ii,jj)=fcars(ii,jj) endif c c read-in dipole derivatives if(lfdp(2))then call READTEN(nmax,P,AAT,NAT,tnamel(i)) c transform dipole derivatives to local system: write(6,*)'dipole der in local system:' do 51 ii=1,N4 ia=i4(i,ii) call trp(nmax,P,e0,ia) do 51 k=1,3 do 51 l=1,3 write(6,643)ii,k,l,P(ia,k,l) 643 format(3i3,f12.3) 51 qqx(i,ii,k,l)=P(ia,k,l) endif c c read-in polarizability derivatives if(lfdp(3))then lcompl=detect(pnamel(i)) if(lcompl)then write(6,*)'complex pol der in local system:' else write(6,*)'pol der in local system:' endif call READTTT(alpha,AT,G,nat,lcompl,pnamel(i)) CALL TTTT(alpha,AT,G,NAT,1,r,lcompl) do 52 ii=1,N4 ia=i4(i,ii) call tra(ia,nat,alpha,qt,e0,lcompl) call tra(ia,nat,G,qg,e0,lcompl) call trt(ia,nat,AT,qa,e0,lcompl) do 52 mm=1,3 do 52 k=1,3 do 52 l=1,3 write(6,644)ii,mm,k,l,qt(mm+3*(ia-1),k,l) 644 format(4i3,f12.3) if(lcompl)then qpx(i,n4+ii,mm,k,l)=qt(3*nat+mm+3*(ia-1),k,l) qgx(i,n4+ii,mm,k,l)=qg(3*nat+mm+3*(ia-1),k,l) ax(i,n4+ii,mm,k,l,1)=qa(3*nat+mm+3*(ia-1),k,l,1) ax(i,n4+ii,mm,k,l,2)=qa(3*nat+mm+3*(ia-1),k,l,2) ax(i,n4+ii,mm,k,l,3)=qa(3*nat+mm+3*(ia-1),k,l,3) endif qpx(i,ii,mm,k,l)=qt(mm+3*(ia-1),k,l) qgx(i,ii,mm,k,l)=qg(mm+3*(ia-1),k,l) ax(i,ii,mm,k,l,1)=qa(mm+3*(ia-1),k,l,1) ax(i,ii,mm,k,l,2)=qa(mm+3*(ia-1),k,l,2) 52 ax(i,ii,mm,k,l,3)=qa(mm+3*(ia-1),k,l,3) endif c 2 deallocate(alpha,AT,G,qt,qg,r,qa,FCAR) c ========================================================= c c set-up matrix (A B = C, B=A^-1 C) nd=n6+1 allocate(a(nd,nd),ipiv(nd),work(3*nd*nd),am(nd,nd),c(nd),d(nd,5), 1b(3*n4,3*n4,nd),bx(n4,3,nd,3),bp(2*n4,3,n6,3,3), 1ba(2*n4,3,n6,3,3,3),bg(2*n4,3,n6,3,3)) a=0.0d0 do 10 i1=1,n6 do 10 j1=1,n6 do 10 i=1,m 10 a(i1,j1)=a(i1,j1)+fi(i,i1)*fi(i,j1) do 12 i1=1,n6 a(i1,nd)=1.0d0 12 a(nd,i1)=1.0d0 c c invert nd0=nd call inv(a,nd0,nd,ipiv,work,am) c b=0.0d0 if(lfdp(1))then do 13 i1=1,3*n4 do 13 j1=1,i1 c=0.0d0 do 14 j=1,n6 do 14 jj=1,m 14 c(j)=c(j)+(ff(jj,i1,j1)-ff(1,i1,j1))*fi(jj,j) do 13 ii=1,n6 do 13 jj=1,nd 13 b(i1,j1,ii)=b(i1,j1,ii)+am(ii,jj)*c(jj) c c check fit write(6,*)' Reference, current and fitted diagonal constants:' do 18 i=1,m do 18 j=1,7 sum=ff(1,j,j) do 19 ii=1,n6 19 sum=sum+fi(i,ii)*b(j,j,ii) 18 write(6,1605)i,j,ff(1,j,j),ff(i,j,j),sum 1605 format(2i3,3f12.6) endif c bx=0.0d0 if(lfdp(2))then c c form c for each eigenvalue and dipole component, dipole FIT do 131 i=1,n4 do 131 iy=1,3 do 131 ix=1,3 c=0.0d0 do 141 j=1,n6 do 141 jj=1,m 141 c(j)=c(j)+(qqx(jj,i,iy,ix)-qqx(1,i,iy,ix))*fi(jj,j) do 131 ii=1,n6 do 131 jj=1,nd 131 bx(i,iy,ii,ix)=bx(i,iy,ii,ix)+am(ii,jj)*c(jj) write(6,*)' Old and fitted dipole constants:' do 181 i=1,m do 181 j=1,4 do 181 jx=1,3 do 181 ix=1,3 old=qqx(i,j,jx,ix)-qqx(1,j,jx,ix) sum=0.0d0 do 191 ii=1,n6 191 sum=sum+fi(i,ii)*bx(j,jx,ii,ix) 181 write(6,16051)i,j,jx,ix,qqx(i,j,jx,ix),old,sum 16051 format(4i3,3f12.5) endif c bp=0.0d0 bg=0.0d0 ba=0.0d0 if(lfdp(3))then c if(lcompl)then nc=0 else nc=1 endif c form c for each component, polarizability FIT do 132 ic=0,nc do 132 i=1,n4 do 132 iy=1,3 do 132 ix=1,3 do 132 iz=1,3 d=0.0d0 do 142 j=1,n6 do 142 jj=1,m d(j,1)=d(j,1)+ 1 (qpx(jj,i+n4*ic,iy,ix,iz)-qpx(1,i+n4*ic,iy,ix,iz))*fi(jj,j) d(j,2)=d(j,2)+ 1 (qgx(jj,i+n4*ic,iy,ix,iz)-qgx(1,i+n4*ic,iy,ix,iz))*fi(jj,j) d(j,3)=d(j,3)+ 1 (ax(jj,i+n4*ic,iy,ix,iz,1)-ax(1,i+n4*ic,iy,ix,iz,1))*fi(jj,j) d(j,4)=d(j,4)+ 1 (ax(jj,i+n4*ic,iy,ix,iz,2)-ax(1,i+n4*ic,iy,ix,iz,2))*fi(jj,j) 142 d(j,5)=d(j,5)+ 1 (ax(jj,i+n4*ic,iy,ix,iz,3)-ax(1,i+n4*ic,iy,ix,iz,3))*fi(jj,j) c c get Bp do 132 ii=1,n6 do 132 jj=1,nd bp(i+ic*n4,iy,ii,ix,iz)=bp(i+ic*n4,iy,ii,ix,iz)+am(ii,jj)*d(jj,1) bg(i+ic*n4,iy,ii,ix,iz)=bg(i+ic*n4,iy,ii,ix,iz)+am(ii,jj)*d(jj,2) ba(i+ic*n4,iy,ii,ix,iz,1)=ba(i+ic*n4,iy,ii,ix,iz,1)+am(ii,jj) 1 *d(jj,3) ba(i+ic*n4,iy,ii,ix,iz,2)=ba(i+ic*n4,iy,ii,ix,iz,2)+am(ii,jj) 1 *d(jj,4) 132 ba(i+ic*n4,iy,ii,ix,iz,3)=ba(i+ic*n4,iy,ii,ix,iz,3)+am(ii,jj) 1 *d(jj,5) write(6,*)' Reference, current, and fitted polarizability' do 182 i=1,m do 182 j=1,4 do 182 jx=1,3 do 182 ix=1,3 do 182 iz=1,3 sum=0.0d0 do 192 ii=1,n6 192 sum=sum+fi(i,ii)*bp(j,jx,ii,ix,iz) 182 write(6,16052)i,j,jx,ix,iz,qpx(1,j,jx,ix,iz),qpx(i,j,jx,ix,iz), 1 qpx(1,j,jx,ix,iz)+sum 16052 format(5i3,3f12.5) endif c c save to disk open(55,file='BIJ.SCR',form='unformatted') write(55)n4,n6 write(55)lfdp write(55)fi_ref if(lfdp(1))write(55)b if(lfdp(2))write(55)bx if(lfdp(3))then write(55)bp write(55)bg write(55)ba endif close(55) call report('Fitting constants have been saved to BIJ.SCR.') 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 function detect(s) integer*4 J,i logical detect character*(*) s character*120 s120 open(15,file=s) do 1 J=1,6 1 READ(15,120)s120 120 format(a120) close(15) i=0 do 2 J=1,len(s120) 2 if(s120(J:J).eq.'.')i=i+1 detect=i.gt.3 return end subroutine rg(n,nat,ity,r,ic) implicit none character*(*) n integer*4 nat,ity(nat),ic,ia,ix real*8 r(nat,3) OPEN(20,FILE=n,STATUS='OLD') read(20,*) read(20,*)nat if(ic.eq.1)then do 4 ia=1,nat 4 read(20,*)ity(ia),(r(ia,ix),ix=1,3) endif close(20) return end subroutine rgn(n,nat) implicit none character*(*) n integer*4 nat OPEN(20,FILE=n,STATUS='OLD') read(20,*) read(20,*)nat close(20) return end subroutine lcs(e0,r,nmax,i4,n4,m,i) implicit none integer*4 nmax,n4,m,i4(m,n4),ix,i real*8 e0(3,3),r(nmax,3) do 50 ix=1,3 e0(1,ix)=r(i4(i,2),ix)-r(i4(i,1),ix) 50 e0(3,ix)=r(i4(i,3),ix)-r(i4(i,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) return end subroutine trp(nmax,P,e0,ia) implicit none integer*4 nmax,ia,k,j real*8 P(nmax,3,3),e0(3,3) real*8,allocatable:: T(:,:,:) allocate(T(nmax,3,3)) c c first index do 511 k=1,3 do 511 j=1,3 511 T(ia,j,k)=e0(j,1)*P(ia,1,k)+e0(j,2)*P(ia,2,k)+e0(j,3)*P(ia,3,k) c second index do 512 k=1,3 do 512 j=1,3 512 P(ia,j,k)=e0(k,1)*T(ia,j,1)+e0(k,2)*T(ia,j,2)+e0(k,3)*T(ia,j,3) return end subroutine tra(ia,nat,a,qt,e0,lcompl) implicit none integer*4 nat,ia,i0,i1,i2,i3,l,ic,nc,no real*8 a(6*nat,3,3),qt(6*nat,3,3),e0(3,3) logical lcompl if(lcompl)then nc=2 else nc=1 endif do 1 ic=1,nc if(ic.eq.1)then no=0 else no=3*nat endif c c first(atomic) index do 611 i1=1,3 i0=i1+3*(ia-1)+no do 611 i2=1,3 do 611 i3=1,3 qt(i0,i2,i3)=0.0d0 do 611 l=1,3 611 qt(i0,i2,i3)=qt(i0,i2,i3)+e0(i1,l)*a(no+l+3*(ia-1),i2,i3) c c second index do 612 i1=1,3 i0=i1+3*(ia-1)+no do 612 i2=1,3 do 612 i3=1,3 a(i0,i2,i3)=0.0d0 do 612 l=1,3 612 a(i0,i2,i3)=a(i0,i2,i3)+e0(i2,l)*qt(i0,l,i3) c c third index do 613 i1=1,3 i0=i1+3*(ia-1)+no do 613 i2=1,3 do 613 i3=1,3 qt(i0,i2,i3)=0.0d0 do 613 l=1,3 613 qt(i0,i2,i3)=qt(i0,i2,i3)+e0(i3,l)*a(i0,i2,l) 1 continue return end subroutine trt(ia,nat,a,qt,e0,lcompl) implicit none integer*4 nat,ia,i0,i1,i2,i3,l,ic,nc,no,i4 real*8 a(6*nat,3,3,3),qt(6*nat,3,3,3),e0(3,3) logical lcompl if(lcompl)then nc=2 else nc=1 endif do 1 ic=1,nc if(ic.eq.1)then no=0 else no=3*nat endif c c first(atomic) index do 611 i1=1,3 i0=i1+3*(ia-1)+no do 611 i2=1,3 do 611 i3=1,3 do 611 i4=1,3 qt(i0,i2,i3,i4)=0.0d0 do 611 l=1,3 611 qt(i0,i2,i3,i4)=qt(i0,i2,i3,i4)+e0(i1,l)*a(no+l+3*(ia-1),i2,i3,i4) c c second index do 612 i1=1,3 i0=i1+3*(ia-1)+no do 612 i2=1,3 do 612 i3=1,3 do 612 i4=1,3 a(i0,i2,i3,i4)=0.0d0 do 612 l=1,3 612 a(i0,i2,i3,i4)=a(i0,i2,i3,i4)+e0(i2,l)*qt(i0,l,i3,i4) c c third index do 613 i1=1,3 i0=i1+3*(ia-1)+no do 613 i2=1,3 do 613 i3=1,3 do 613 i4=1,3 qt(i0,i2,i3,i4)=0.0d0 do 613 l=1,3 613 qt(i0,i2,i3,i4)=qt(i0,i2,i3,i4)+e0(i3,l)*a(i0,i2,l,i4) c c fourth index do 614 i1=1,3 i0=i1+3*(ia-1)+no do 614 i2=1,3 do 614 i3=1,3 do 614 i4=1,3 a(i0,i2,i3,i4)=0.0d0 do 614 l=1,3 614 a(i0,i2,i3,i4)=a(i0,i2,i3,i4)+e0(i3,l)*qt(i0,i2,i3,l) 1 continue qt=a return end subroutine trf(nat,f,e0,ia) implicit none integer*4 nat,ia,j,i,i3,ii,j3,ja,jj real*8 f(3*nat,3*nat),e0(3,3) real*8,allocatable:: T(:,:) allocate(T(3*nat,3*nat)) c c first index T=0.0d0 do 1 ia=1,nat i3=3*(ia-1) do 1 i=1,3 ii=i+i3 do 1 jj=1,3*nat 1 T(ii,jj)=e0(i,1)*f(1+i3,jj)+e0(i,2)*f(2+i3,jj)+e0(i,3)*f(3+i3,jj) do 2 ii=1,3*nat do 2 ja=1,nat do 2 j=1,3 j3=3*(ja-1) jj=j+j3 2 f(ii,jj)=e0(j,1)*T(ii,1+j3)+e0(j,2)*T(ii,2+j3)+e0(j,3)*T(ii,3+j3) return end FUNCTION EPS(I,J,K) REAL*8 EPS INTEGER*4 I,J,K EPS=0.0d0 IF (I.EQ.1.AND.J.EQ.2.AND.K.EQ.3)EPS= 1.0d0 IF (I.EQ.1.AND.J.EQ.3.AND.K.EQ.2)EPS=-1.0d0 IF (I.EQ.2.AND.J.EQ.3.AND.K.EQ.1)EPS= 1.0d0 IF (I.EQ.2.AND.J.EQ.1.AND.K.EQ.3)EPS=-1.0d0 IF (I.EQ.3.AND.J.EQ.1.AND.K.EQ.2)EPS= 1.0d0 IF (I.EQ.3.AND.J.EQ.2.AND.K.EQ.1)EPS=-1.0d0 RETURN END c ============================================================ SUBROUTINE TTTT(ALPHA,A,G,NAT,ICO,R,lcomplex) c Tensors derivatives C Transforms A and G to local origins (ICO=1) or c common origin (ICO=2); C sets A and G to zero for ICO=3. C R supposed to be passed in in A C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION ALPHA(6*NAT,3,3),G(6*NAT,3,3),A(6*NAT,3,3,3), 1R(NAT,3) real*8,allocatable::X(:,:) logical lcomplex IF(ICO.EQ.3)then G=0.0d0 A=0.0d0 WRITE(3,*)' A,G set to zero' return ENDIF if(lcomplex)then M=1 else M=0 endif C C transfer into atomic coordinates: allocate(X(3,NAT)) DO 7 I=1,3 DO 7 J=1,NAT 7 X(I,J)=R(J,I)/0.529177D0 C IF(ICO.EQ.1)SIGN=1.0d0 IF(ICO.EQ.2)SIGN=-1.0d0 W=1.0d0 C supposedly the static limit G/w is calculated C DO 2 IA=1,3 DO 2 IB=1,3 DO 2 L=1,NAT DO 2 IE=1,3 if(lcomplex)then i1=3*(L-1)+IE i2=3*(L-1)+IE+3*NAT SUM=0.0d0 DO 5 IG=1,3 DO 5 ID=1,3 5 SUM=SUM+SIGN*(0.5d0*W)*EPS(IB,IG,ID)*X(IG,L)*ALPHA(i1,IA,ID) G(i2,IA,IB)=G(i2,IA,IB)+SUM SUM=0.0d0 DO 6 IG=1,3 DO 6 ID=1,3 6 SUM=SUM+SIGN*(0.5d0*W)*EPS(IB,IG,ID)*X(IG,L)*ALPHA(i2,IA,ID) G(i1,IA,IB)=G(i1,IA,IB)+SUM else I1=3*(L-1)+IE SUM=0.0d0 DO 1 IG=1,3 DO 1 ID=1,3 1 SUM=SUM+SIGN*(0.5d0*W)*EPS(IB,IG,ID)*X(IG,L)*ALPHA(i1,IA,ID) G(i1,IA,IB)=G(i1,IA,IB)+SUM endif 2 CONTINUE C c A IA index is electric DO 3 IA=1,3 DO 3 IB=1,3 DO 3 IC=1,3 DO 3 L=1,NAT DO 3 IE=1,3 DO 3 icc=0,M IIND=3*(L-1)+IE+3*NAT*icc SUM=0.0d0 IF(IB.EQ.IC)THEN DO 4 ID=1,3 4 SUM=SUM+X(ID,L)*ALPHA(IIND,IA,ID)*SIGN ENDIF A(IIND,IA,IB,IC)=A(IIND,IA,IB,IC)+SUM 1-1.5d0*SIGN*(X(IB,L)*ALPHA(IIND,IA,IC)+X(IC,L)*ALPHA(IIND,IA,IB)) 3 CONTINUE IF(ICO.EQ.1)WRITE(3,*)' A,G transformed into atomic origins' IF(ICO.EQ.2)WRITE(3,*)' A,G transformed into common origin' RETURN END SUBROUTINE READTTT(ALPHA,A,G,NAT,lcomplex,s) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION ALPHA(6*NAT,3,3),G(6*NAT,3,3),A(6*NAT,3,3,3) logical lcomplex character*(*) s open(15,file=s,status='old') M=0 if(lcomplex)M=3 READ(15,*) READ(15,*)NATC if(NAT.NE.NATC)call report('Wrong number of atoms in .TTT file') READ(15,*) READ(15,*) DO 1 I=1,3 READ(15,*) DO 1 L=1,NAT DO 1 IX=1,3 IIND=3*(L-1)+IX 1 READ(15,*)(ALPHA(IIND,I,J),J=1,2),(ALPHA(IIND,I,J),J=1,3), 1(ALPHA(3*NAT+IIND,I,J),J=1,M) READ(15,*) READ(15,*) DO 2 I=1,3 READ(15,*) DO 2 L=1,NAT DO 2 IX=1,3 IIND=3*(L-1)+IX 2 READ(15,*)(G(IIND,I,J),J=1,2),(G(IIND,I,J),J=1,3), 1(G(3*NAT+IIND,I,J),J=1,M) READ(15,*) READ(15,*) DO 3 I=1,3 DO 3 J=1,3 READ(15,*) DO 3 L=1,NAT DO 3 IX=1,3 IIND=3*(L-1)+IX 3 READ(15,*)(A(IIND,I,J,K),K=1,2),(A(IIND,I,J,K),K=1,3), 1(A(3*NAT+IIND,I,J,K),K=1,M) CLOSE(15) RETURN END c SUBROUTINE READTEN(N0,P,A,NAT,s) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION P(N0,3,3),A(N0,3,3) character*(*) s OPEN(15,FILE=s,STATUS='OLD') 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,*) close(15) RETURN END C SUBROUTINE READFF(N,FCAR,s) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(N,N) character*(*) s OPEN(20,FILE=s,STATUS='OLD') 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 DO 31 I=1,N DO 31 J=I+1,N 31 FCAR(I,J)=FCAR(J,I) WRITE(*,*)' Cartesian FF read in ... ' CLOSE(20) RETURN END C c ====================================================================== subroutine inv(fx,M,M0,ipiv,work,a) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION fx(M0,M0),a(M0,M0),ipiv(M0),work(3*M0*M0) do 3 i=1,M do 3 j=i,M a(j,i)=fx(i,j) 3 a(i,j)=fx(i,j) M1=M M2=M call dgetrf(M1,M2,a,M0,ipiv,info) if(info.eq.0)then write(6,*)' Pivoting OK' else call report('pivoting failed') endif c lwork=3*M*M call dgetri(M,a,M0,ipiv,work,lwork,info) c if(info.eq.0)write(6,*)'successful inversion' if(info.lt.0)then