program getqd c get quadrupole derivatives from FILE.OUT implicit none integer*4 nat,i,ir,iu,iq,ia,ix,id,j,ip,ii integer*4 ,allocatable::uu(:) real*8 u(3),u0(3),q(6),q0(6),tol,step,du,dq,dp real*8,allocatable::r(:),r0(:),ud(:),qd(:) character*1 xyz(3) data xyz/'x','y','z'/ nat=0 tol=0.0001d0 open(9,file='FILE.OUT',status='old') c read number of atoms: call readf(nat,1,r0,u0,q0,ir,iu,iq) allocate(r(3*nat),r0(3*nat),ud(9*nat),qd(18*nat),uu(3*nat)) c read first geometry and moments: call readf(nat,2,r0,u0,q0,ir,iu,iq) if(ir.ne.0)then write(6,*)'Zero point, coordinates read' else call report('R not found') endif if(iu.ne.0)then write(6,*)'Zero point dipole X Y Z:' write(6,601)(u0(ix),ix=1,3) 601 format(6f12.6) else call report('U not found') endif if(iq.ne.0)then write(6,*)'Zero point quadrupole XX YY ZZ XY XZ YZ:' write(6,601)(q0(ix),ix=1,6) else call report('Q not found') endif do 3 i=1,9*nat if(i.le.3*nat)uu(i)=0 ud(i)=0.0d0 qd(i )=0.0d0 3 qd(i+9*nat)=0.0d0 write(6,*)6*nat,' points expected' do 1 ip=1,6*nat call readf(nat,2,r,u,q,ir,iu,iq) if(ir.eq.0)call report('R not found') if(iu.eq.0)call report('U not found') if(iq.eq.0)call report('Q not found') ia=0 ix=0 id=0 do 2 i=1,nat do 2 j=1,3 if(dabs(r(j+3*(i-1))-r0(j+3*(i-1))).gt.tol)then ia=i ix=j id=id+1 endif 2 continue if(id.eq.0)call report('no differentiation') if(id.gt.1)call report('unclear point') ii=ix+3*(ia-1) step=r(ii)-r0(ii) write(6,602)ip,ia,xyz(ix),step 602 format(' Point ',i6,', atom ',i4,' ',a1,f10.5) uu(ii)=uu(ii)+1 do 4 i=1,3 du=(u(i )-u0(i ))/step dq=(q(i )-q0(i ))/step dp=(q(i+3)-q0(i+3))/step if(uu(ii).eq.1)then ud(ii+3*nat*(i-1))=du qd(ii+3*nat*(i-1))=dq qd(ii+3*nat*(i+2))=dp else if(uu(ii).gt.2)call report('IU > 2') ud(ii+3*nat*(i-1))=(ud(ii+3*nat*(i-1))+du)/2.0d0 qd(ii+3*nat*(i-1))=(qd(ii+3*nat*(i-1))+dq)/2.0d0 qd(ii+3*nat*(i+2))=(qd(ii+3*nat*(i+2))+dp)/2.0d0 endif 4 continue 1 continue close(9) write(6,63)(uu(i),i=1,3*nat) 63 format(50i1) call sumr(qd,nat,xyz,u0) c call wrqd('FILE.RAW.QEN',q0,qd,nat,xyz) c call project(qd,nat,r0) c call sumr(qd,nat,xyz,u0) call wrqd('FILE.QEN',q0,qd,nat,xyz) call wrud('FILE.NUM.TEN',ud,nat) end SUBROUTINE wrud(fn,ud,nat) IMPLICIT none character*(*) fn integer*4 nat,i,ia,ix,ic REAL*8 ud(*),debye,bohr REAL*8,allocatable::d(:) allocate(d(9*nat)) debye=2.541746473d0 bohr=0.529177d0 do 2 i=1,9*nat 2 d(i)=ud(i)/debye*bohr OPEN(15,FILE=fn) WRITE(15,1500) NAT,NAT-6,0 1500 FORMAT(3I5) DO 10 ia=1,nat DO 10 ix=1,3 10 WRITE(15,1501) (d(ix+3*(ia-1)+3*nat*(i-1)),i=1,3),ia 1501 FORMAT(3F14.8,I5) DO 220 ic=1,3 DO 220 ia=1,NAT DO 220 ix=1,3 220 WRITE(15,1501) (0.0d0,i=1,3),ia CLOSE(15) WRITE(6,*)fn//' written' RETURN END C subroutine sumr(qd,nat,xyz,u) implicit none character*1 xyz(*) real*8 qd(*),s(6),u(*) integer*4 i,ia,ix,nat do 1 ix=1,3 do 2 i=1,6 s(i)=0.0d0 do 2 ia=2,nat 2 s(i)=s(i)-qd(ix+3*(ia-1)+3*nat*(i-1)) c XX YY ZZ XY XZ YZ c 1 2 3 4 5 6 if(ix.eq.1)then s(1)=s(1)+u(1)+u(1) s(4)=s(4)+u(2) s(5)=s(5)+u(3) endif if(ix.eq.2)then s(2)=s(2)+u(2)+u(2) s(4)=s(4)+u(1) s(6)=s(6)+u(3) endif if(ix.eq.3)then s(3)=s(3)+u(3)+u(3) s(5)=s(5)+u(1) s(6)=s(6)+u(2) endif 1 write(6,60)xyz(ix),(qd(ix+3*nat*(i-1)),i=1,6),s 60 format(' Sum rules ',a1,/,6f12.4,/,6f12.4,/) return end subroutine wrqd(fn,q0,qd,nat,xyz) implicit none character*(*) fn character*1 xyz(*) integer*4 i,ia,ix,nat real*8 q0(*),qd(*),bohr,debye real*8,allocatable::q(:),d(:) bohr=0.529177d0 debye=2.541746473d0 allocate(q(6),d(18*nat)) do 1 i=1,6 1 q(i)=q0(i)/debye/bohr do 3 i=1,18*nat 3 d(i)=qd(i)/debye open(9,file=fn) write(9,90)nat,'XX',q(1),'YY',q(2),'ZZ',q(3), 1 'XY',q(4),'XZ',q(5),'YZ',q(6) 90 format(i6,' atoms',/, 1 ' Equilibrium quadrupole, sum(i)qi ri ri, atomic units:',/, 12(3(3x,a2,2x,f19.4),/)) write(9,91) 91 format(' Quadrupole derivatives, atomic units: ',/, 1 ' Atom Coord XX YY ZZ ',/, 1 ' XY XZ YZ ',/, 1 ' ------------------------------------------------') do 2 ia=1,nat do 2 ix=1,3 if(ix.eq.1)then write(9,92)ia,xyz(ix),(d(ix+3*(ia-1)+3*nat*(i-1)),i=1,6) 92 format(i5,3x,a1,3e13.4,/,9x,3e13.4) else write(9,93)xyz(ix),(d(ix+3*(ia-1)+3*nat*(i-1)),i=1,6) 93 format(8x,a1,3e13.4,/,9x,3e13.4) endif 2 continue close(9) write(6,*)fn//' written' return end subroutine readf(nat,ic,r,u,q,ir,iu,iq) implicit none integer*4 nat,i,ic,ix,ir,iu,iq real*8 r(*),u(*),q(*) character*80 s80 ir=0 iu=0 iq=0 c c read number of atoms and rewind: if(ic.eq.1)then 1 read(9,90,end=99,err=99)s80 90 format(a80) if(s80(27:44).eq.'Input orientation:'.and.nat.eq.0)then do 2 i=1,4 2 read(9,*) 3 read(9,90,end=88,err=88)s80 if(s80(2:3).ne.'--')then nat=nat+1 goto 3 endif write(6,*)nat,' atoms' rewind 9 return endif goto 1 99 close(9) return endif c c read everything if(ic.eq.2)then 4 read(9,90,end=88,err=88)s80 if(s80(27:44).eq.'Input orientation:')then do 5 i=1,4 5 read(9,*) do 6 i=1,nat 6 read(9,*)(r(ix+3*(i-1)),ix=1,3),(r(ix+3*(i-1)),ix=1,3) ir=1 endif if(s80(2:39).eq.'Dipole moment (field-independent basis')then read(9,91)(u(ix),ix=1,3) 91 format(3(7x,f19.4)) iu=1 endif if(s80(2:39).eq.'Quadrupole moment (field-independent b')then read(9,91)(q(ix),ix=1,3) read(9,91)(q(ix),ix=4,6) c order: XX YY ZZ XY XZ YZ iq=1 return endif goto 4 endif 88 close(9) call report('End of file') end subroutine report(s) character*(*) s write(6,*)s stop end subroutine project(qd,nat,r0) implicit none integer*4 nat,N,N6,i,ia,ix,j real*8 qd(*),r0(*) real*8,allocatable::R(:,:),A(:,:),t(:) N=3*nat N6=6 allocate(R(3,NAT),A(N,N),t(N)) do 3 ia=1,nat do 3 ix=1,3 3 R(ix,ia)=r0(ix+3*(ia-1)) call DOMA(A,R,NAT,N6) do 1 ix=1,6 do 2 i=1,N t(i)=0.0d0 do 2 j=1,N 2 t(i)=t(i)+A(i,j)*qd(j+N*(ix-1)) do 1 i=1,N 1 qd(i+N*(ix-1))=t(i) write(6,*)' translation/rotation projected' return end SUBROUTINE DOMA(A,C,NAT,N6) c N6 - number of coordinates to be projected c 1..3 translations, 4..6 rotations IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION C(3,NAT),A(3*NAT,3*NAT) real*8,allocatable::TEM(:,:) allocate(TEM(3*NAT,3*NAT)) N=NAT*3 AMACH=1.0d-12 DO 12 I=1,N DO 12 J=1,N 12 A(I,J)=0.0d0 DO 1 I=1,NAT A(3*I-2,1)= 1.00d0 A(3*I-1,2)= 1.00d0 A(3*I ,3)= 1.00d0 A(3*I-2,4)=-C(2,I) A(3*I-1,4)= C(1,I) A(3*I-1,5)=-C(3,I) A(3*I ,5)= C(2,I) A(3*I-2,6)= C(3,I) 1 A(3*I ,6)=-C(1,I) CALL ORT(A,N,N6,NAT) ICOL=1 5 IF(ICOL.GT.N6)GOTO 8 S=0.0d0 DO 6 J=1,N 6 S=S+A(J,ICOL)**2 IF(S.LT.1000.0d0*AMACH)THEN DO 7 J=1,N 7 A(J,ICOL)=A(J,N6) N6=N6-1 GOTO 5 ENDIF ICOL=ICOL+1 GOTO 5 8 CONTINUE WRITE(*,*)N6,' coordinates projected' DO 2 I=1,N DO 2 J=1,N TEM(I,J)=0.0d0 DO 2 II=1,N6 2 TEM(I,J)=TEM(I,J)+A(I,II)*A(J,II) DO 3 I=1,N DO 4 J=1,N 4 A(I,J)=-TEM(I,J) 3 A(I,I)=A(I,I)+1.0d0 RETURN END c =============================== SUBROUTINE ORT(A,NDIM,N6,NAT) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION A(NDIM,NDIM) AMACH=1.0d-13 NAT3=3*NAT C=0.0d0 DO 1 I=1,NAT3 1 C=C+A(I,1)**2 IF(C.GT.AMACH)C=1.0d0/SQRT(C) DO 2 I=1,NAT3 2 A(I,1)=A(I,1)*C DO 3 IC=2,N6 C Column IC, orthonormalize to previous columns: DO 3 ICP=1,IC-1 SP=0.0d0 DO 5 J=1,NAT3 5 SP=SP+A(J,IC)*A(J,ICP) C Subtract projection to ICP: DO 6 J=1,NAT3 6 A(J,IC)=A(J,IC)-SP*A(J,ICP) C Normalize : SP=0.0d0 DO 7 J=1,NAT3 7 SP=SP+A(J,IC)**2 IF(SP.GT.AMACH)SP=1.0d0/SQRT(SP) DO 3 J=1,NAT3 3 A(J,IC)=A(J,IC)*SP RETURN END