MODULE tmod TYPE :: itensors real*8,allocatable:: alphad(:,:,:),gpd(:,:,:),ad(:,:,:,:) real*8,allocatable:: tw(:) END TYPE itensors contains SUBROUTINE READTTTQ(ft,nm,ite,ie) IMPLICIT none INTEGER*4 nm,nq,LL,I,J,K,ie character*80 ft TYPE(itensors)::ite(*) open(15,file=ft) READ(15,*) READ(15,*)nq if(nq.ne.nm)call report('nq<>ntw') READ(15,*) READ(15,*) DO 1 LL=1,nq READ(15,*)ite(ie)%tw(LL),ite(ie)%tw(LL) DO 1 I=1,3 1 READ(15,*)ite(ie)%alphad(LL,I,1),(ite(ie)%alphad(LL,I,J),J=1,3) ite(ie)%tw(nq+1)=0.0d0 READ(15,*) READ(15,*) DO 2 LL=1,abs(nq) READ(15,*) DO 2 I=1,3 2 READ(15,*)ite(ie)%GPD(LL,I,1),(ite(ie)%GPD(LL,I,J),J=1,3) READ(15,*) READ(15,*) DO 3 LL=1,abs(nq) READ(15,*) DO 3 I=1,3 DO 3 J=1,3 3 READ(15,*)(ite(ie)%Ad(LL,I,J,K),K=1,2), 1(ite(ie)%Ad(LL,I,J,K),K=1,3) c first index electric CLOSE(15) RETURN END SUBROUTINE TTT1(alpha,A,G,ICO,r,n,ite,lusea, 1luseg,ntw,ltr) c Tensors, not derivatives C Transforms A and G to local origins (ICO=1), c common origin (ICO=2), C sets them to zero (ICO=3). C c supposed to be passed in aus C the static limit G/w is in G C implicit none INTEGER*4 ICO,n,iu,IA,IB,IG,ID,IC,ntw(*),iw real*8 alpha(n,3,3),G(n,3,3),A(n,3,3,3),SIGN,SUM,EPS, 1r(n,3) logical luseg,lusea,ltr(*) TYPE(itensors)::ite(*) c IF(ICO.EQ.1)SIGN= 1.0d0 IF(ICO.EQ.2)SIGN=-1.0d0 IF(ICO.EQ.3)SIGN= 0.0d0 C do 99 iu=1,n if(ltr(iu))then C DO 2 IA=1,3 DO 2 IB=1,3 SUM=0.0d0 DO 1 IG=1,3 DO 1 ID=1,3 1 SUM=SUM+SIGN*0.5d0*EPS(IB,IG,ID)*r(iu,IG)*alpha(iu,IA,ID) G(iu,IA,IB)=G(iu,IA,IB)+SUM IF(ICO.EQ.3.or..not.luseg)G(iu,IA,IB)=0.0d0 do 21 iw=1,ntw(iu) sum=0.0d0 DO 11 IG=1,3 DO 11 ID=1,3 11 SUM=SUM+SIGN*0.5d0*EPS(IB,IG,ID)*r(iu,IG)* 1ite(iu)%alphad(iw,IA,ID) ite(iu)%Gpd(iw,IA,IB)=ite(iu)%Gpd(iw,IA,IB)+SUM c G last index (IB) magnetic 21 IF(ICO.EQ.3.or..not.luseg)ite(iu)%Gpd(iw,IA,IB)=0.0d0 2 CONTINUE C c A .. first index (IA) is electric DO 3 IA=1,3 DO 3 IB=1,3 DO 3 IC=1,3 SUM=0.0d0 IF(IB.EQ.IC)SUM=r(iu,1)*alpha(iu,IA,1)+ 1r(iu,2)*alpha(iu,IA,2)+r(iu,3)*alpha(iu,IA,3) A(iu,IA,IC,IB)=A(iu,IA,IC,IB)+SIGN* 1(SUM-1.5d0*(r(iu,IB)*alpha(iu,IA,IC)+r(iu,IC)*alpha(iu,IA,IB))) IF(ICO.EQ.3.or..not.lusea)A(iu,IA,IC,IB)=0.0d0 do 31 iw=1,ntw(iu) SUM=0.0d0 IF(IB.EQ.IC)SUM= ite(iu)%alphad(iw,IA,1)*r(iu,1) 1 + ite(iu)%alphad(iw,IA,2)*r(iu,2) 2 + ite(iu)%alphad(iw,IA,3)*r(iu,3) ite(iu)%Ad(iw,IA,IC,IB)=ite(iu)%Ad(iw,IA,IC,IB)+ 1SIGN*(SUM-1.5d0*( ite(iu)%alphad(iw,IA,IC)*r(iu,IB) 2 + ite(iu)%alphad(iw,IA,IB)*r(iu,IC))) 31 IF(ICO.EQ.3.or..not.lusea)ite(iu)%Ad(iw,IA,IC,IB)=0.0d0 3 CONTINUE write(6,600)iu 600 format(' molecule ',i4,' transformed') else write(6,601)iu 601 format(' molecule ',i4,' not transformed') endif 99 continue RETURN END function res(iy,ite,it,i) implicit none real*8 res integer*4 iy,it,i,ix,jx,kx,iyo TYPE(itensors)::ite(*) iyo=iy kx=(iy-1)/1000 iy=iy-1000*kx jx=(iy-1)/100 iy=iy-100*jx ix=(iy-1)/10 iy=iy-10*ix if(iy.eq.1)then res=ite(i)%alphad(it,ix,jx) else if(iy.eq.2)then res=ite(i)%gpd(it,ix,jx) else if(iy.eq.4)then res=-ite(i)%gpd(it,ix,jx) else if(iy.eq.3)then res=ite(i)%ad(it,ix,jx,kx) else write(6,*)iyo,iy,ix,jx,kx call report('error in res') endif endif endif endif return end subroutine MDE(n,ndim,x,Itr,Pt,ntw,tolw,ite,w) implicit none integer*4 ndim,i,j,ii,it,jj,jt,kk,n,ntw(*),iy, 1Itr(ndim,ndim),ko,ll real*8 x(n*24,n*24),w,tolw,wj, 1Pt(ndim,ndim),sum,wi real*8 tem(24) TYPE(itensors)::ite(*) c make X.P, associated transition frequency list is with j-index: c make total polarizability Pt = P(E-X.P)^(-1) for W <- 0: c = Ptr_w w' temi_w'_0= P . (1+XP)= c P.tem: do 74 j=1,n do 74 jt=1,ntw(j) wj=ite(j)%tw(jt) do 74 i=1,n do 74 it=1,ntw(i) wi=ite(i)%tw(it) if(dabs(w-wi-wj).lt.tolw)then do 76 jj=24*(j-1)+1,24*j c prepare element (E+XP)_jt,kk,jj do 1 kk=1,24 ko=24*(i-1)+kk sum=0.0d0 do 2 ll=24*(j-1)+1,24*j iy=Itr(ll,jj) 2 if(iy.ne.0)sum=sum+x(ko,ll)*res(iy,ite,jt,j) if(jt.eq.ntw(j).and.ko.eq.jj)sum=sum+1.0d0 1 tem(kk)=sum c perturbation approach, (E-X.P)^-1 ~ (E+X.P) do 76 ii=24*(i-1)+1,24*i c P.(E+X.P): sum=0.0d0 do 761 kk=1,24 ko=24*(i-1)+kk c actually kk=1,ndim, but Ptr diagonal in i iy=Itr(ii,ko) 761 if(iy.ne.0)sum=sum+res(iy,ite,it,i)*tem(kk) 76 Pt(ii,jj)=Pt(ii,jj)+sum endif 74 continue return end END MODULE tmod c ============================================================ PROGRAM IRROA use tmod IMPLICIT none integer*4 n,i,ix,jx,kx,ii,istart,nm,NW,k,ntt,ic,j, 1jt,kt,ktend,l,lt,ltend,lend,is1,is2,isu,is2s, 1it,ndim,iw integer*4 ,allocatable::ntw(:),Itr(:,:) real*8 alphar(3,3),gpr(3,3),ar(3,3,3) real*8,allocatable:: alpha(:,:,:),gp(:,:,:),a(:,:,:,:), 1r(:,:),X(:,:),rams(:),roas(:), 1wl(:),Pt(:,:) real*8 tol, 1WMIN,WMAX,DW,GAMMA,W0,TEMP,tempcm,w,tolw character*80,allocatable:: fn(:),ftw(:) logical,allocatable:: ltr(:) character*80 ft,s80 logical lave,luseg,lusea,lcont,lff,lgauss TYPE(itensors),allocatable::ite(:) c lcont ... try the "continuum" frequency model c lff .. try also frequency combinations wi+wj lff=.false. lgauss=.false. luseg=.true. lusea=.true. lcont=.false. tol=0.01d0 lave=.false. WMAX=3500.0d0 WMIN=0.0d0 W0=0.0d0 NW=1751 GAMMA=10.0d0 TEMP=273.0d0 tolw=0.2d0 c iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii c read keywords: open(9,file='IRROA.OPT',status='old') 66 read(9,920,end=888,err=888)s80 920 format(a80) if(s80(1:4).eq.'AVER')read(9,*)lave if(s80(1:4).eq.'ADER')read(9,*)lusea if(s80(1:4).eq.'GDER')read(9,*)luseg if(s80(1:3).eq.'TOL')read(9,*)tol if(s80(1:4).eq.'WMAX')read(9,*)WMAX if(s80(1:4).eq.'WTOL')read(9,*)tolw if(s80(1:4).eq.'WMIN')read(9,*)WMIN if(s80(1:4).eq.'TEMP')read(9,*)TEMP if(s80(1:4).eq.'TOLW')read(9,*)tolw if(s80(1:3).eq.'LFF')read(9,*)lff if(s80(1:4).eq.'LCON')read(9,*)lcont if(s80(1:4).eq.'LGAU')read(9,*)lgauss if(s80(1:2).eq.'W0')read(9,*)W0 if(s80(1:2).eq.'NW')read(9,*)NW if(s80(1:5).eq.'GAMMA')read(9,*)GAMMA goto 66 888 close(9) c c read again for filnames: open(9,file='IRROA.OPT',status='old') c number of particles/chromophores: read(9,*)n allocate(r(n,3),fn(n),ftw(n),ntw(n),ltr(n), 1alpha(n,3,3),gp(n,3,3),a(n,3,3,3),ite(n)) do 1 i=1,n ltr(i)=.true. read(9,900)fn(i) 900 format(a80) c geometries, get the centers only: 1 CALL READRX(r,n,fn,i) c normal mode polarizability derivatives c 1) read names and find maximum number of modes: nm=0 write(6,*)' group number of transition frequencies:' do 2 i=1,n read(9,900)ftw(i) open(15,file=ftw(i)) READ(15,*) c ntw: number of transition frequencies: READ(15,*)ntw(i) write(6,*)i,ntw(i) close(15) allocate(ite(i)%alphad(ntw(i)+1,3,3),ite(i)%gpd(ntw(i)+1,3,3), 1ite(i)%ad(ntw(i)+1,3,3,3),ite(i)%tw(ntw(i)+1)) 2 if(abs(ntw(i)).gt.nm)nm=abs(ntw(i)) c 2) read derivatives/transition polarizabilities/ do 21 i=1,n 21 CALL READTTTQ(ftw(i),ntw(i),ite,i) write(6,*)'1' c polarizability alpha do 3 i=1,n read(9,900)ft 3 call readpol(alpha,i,ft,n) c G': do 4 i=1,n read(9,900)ft 4 call readpolg(gp,i,ft,n) c A: do 5 i=1,n read(9,900)ft 5 call readpola(a,i,ft,n) write(6,*)'1' 67 read(9,920,end=889,err=889)s80 c transform the tensors to the molecular positions: if(s80(1:4).eq.'TRAN')read(9,*)(ltr(i),i=1,n) if(s80(1:4).eq.'SAME')then read(9,*)is1,is2 if(is2.lt.0)then c transition frequencies of groups is1...is2 are the same: is2s=is1+1 else c transition frequencies of group is2 are the same as those of is1: is2s=is2 endif do 68 isu=is2s,abs(is2) if(ntw(is1).ne.ntw(isu))call report('ntw(is1) <> ntw(is2)') do 68 i=1,ntw(isu) 68 ite(isu)%tw(i)=ite(is1)%tw(i) endif goto 67 889 close(9) write(6,*)'tensors read' DW=(WMAX-WMIN)/dble(NW-1) tempcm=0.6950d0*TEMP ndim=24*n write(6,6000)lave,lusea,luseg,tol, 1WMAX,WMIN,TEMP,W0,NW,DW,n,ndim,tolw,lff,lcont,GAMMA 6000 format(' Averaging :',l10,/, 1 ' Use A :',l10,/, 2 ' Use G :',l10,/, 2 ' tol :',G20.6,/, 2 ' Wmax :',F20.2,/, 2 ' Wmin :',F20.2,/, 2 ' Temperature :',F20.2,/, 2 ' Wexc(relative) :',F20.2,/, 2 ' Nw :',I20,/, 2 ' DW :',F20.2,/, 2 ' Number of mols. :',i20,/, 2 ' Total dimension :',i20,/, 2 ' tolw :',f20.4,/, 2 ' lff :',l20,/, 2 ' continuum model :',l20,/, 2 ' Gamma :',F20.2,/) c transform G, A to molecular centers call TTT1(alpha,a,gp,1,r,n,ite,lusea,luseg, 1ntw,ltr) write(6,*)' Tensors transformed to molecular centers' write(6,*)' ndim ',ndim allocate(X(ndim,ndim),Itr(ndim,ndim)) call tz(x,ndim,ndim) do 27 ix=1,ndim do 27 jx=1,ndim 27 Itr(ix,jx)=0 write(6,*)' X,P allocated' c ************* Calculate general distance tensor: ********* call doX(n,r,X) write(6,*)' Distance tensor done' c Itr - index for mapping the big matrix on polarizabilities: c a 0 0 G A/3 0 c 0 a -G 0 0 A/3 c 0 -Gt 0 0 0 0 c Gt 0 0 0 0 0 c A/3 0 0 0 0 0 c 0 A/3 0 0 0 0 do 71 i=1,n istart=24*(i-1)+1 do 71 ix=1,3 ii=0 do 71 jx=1,3 Itr(istart+ix-1,istart+jx-1)= 1+10*ix+100*jx Itr(istart+ix+2,istart+jx+2)= 1+10*ix+100*jx Itr(istart+ix-1,istart+jx+8)= 2+10*ix+100*jx Itr(istart+jx+8,istart+ix-1)= 2+10*ix+100*jx Itr(istart+ix+2,istart+jx+5)= 4+10*ix+100*jx Itr(istart+jx+5,istart+ix+2)= 4+10*ix+100*jx do 71 kx=jx,3 ii=ii+1 Itr(istart+ix- 1,istart+ii+11)=3+10*ix+100*jx+1000*kx Itr(istart+ii+11,istart+ix- 1)=3+10*ix+100*jx+1000*kx Itr(istart+ix+ 2,istart+ii+17)=3+10*ix+100*jx+1000*kx 71 Itr(istart+ii+17,istart+ix+ 2)=3+10*ix+100*jx+1000*kx c ************* Calculate general polarizability tensor: ********* c (ordinary polarizabilities as ntw+1^th component) do 72 i=1,n ntw(i)=ntw(i)+1 do 72 ix=1,3 do 72 jx=1,3 ite(i)%alphad(ntw(i),ix,jx)=alpha(i,ix,jx) ite(i)%gpd(ntw(i),ix,jx)=gp(i,ix,jx) do 72 kx=1,3 do 721 it=1,ntw(i)-1 721 ite(i)%ad(ntw(i),ix,jx,kx)=ite(i)%ad(ntw(i),ix,jx,kx)/3.0d0 72 ite(i)%ad(ntw(i),ix,jx,kx)=a(i,ix,jx,kx)/3.0d0 write(6,*)' Zero frequency polarizability added' c c how many total t-frequencies wi,wi+wj: ntt=0 do 22 ic=1,2 c ic=1 ... only count c ic=2 ... add to list if(ic.eq.2)then allocate(wl(ntt)) ntt=0 endif c look at wi do 17 i=1,n do 17 it=1,ntw(i) w=ite(i)%tw(it) if(w.gt.WMAX)goto 17 c c was it already among wk: do 18 k=1,i ktend=ntw(k) if(k.eq.i)ktend=it-1 do 18 kt=1,ktend 18 if(dabs(w-ite(k)%tw(kt)).lt.tolw)goto 17 c w is unique: ntt=ntt+1 if(ic.eq.2)wl(ntt)=w 17 continue c look at wi+wj: if(lff)then do 19 i=1,n do 19 it=1,ntw(i) do 19 j=1,i-1 do 19 jt=1,ntw(j) w=ite(i)%tw(it)+ite(j)%tw(jt) if(w.gt.WMAX)goto 19 c c was it already among wk: do 20 k=1,n do 20 kt=1,ntw(k) 20 if(dabs(w-ite(k)%tw(kt)).lt.tolw)goto 19 c was it already among wk+wl: do 26 k=1,i ktend=ntw(k) if(k.eq.i)ktend=it do 26 kt=1,ktend lend=k-1 if(k.eq.i.and.kt.eq.it)lend=j do 26 l=1,lend ltend=ntw(l) if(k.eq.i.and.kt.eq.it.and.l.eq.j)ltend=jt-1 do 26 lt=1,ltend 26 if(dabs(w-ite(k)%tw(kt)-ite(l)%tw(lt)).lt.tolw)goto 19 c w is unique ntt=ntt+1 if(ic.eq.2)wl(ntt)=w 19 continue endif 22 continue write(6,6012)ntt 6012 format(i9,' unique frequencies:') write(6,6013)(wl(i),i=1,ntt) 6013 format(8F10.2) write(6,*) allocate(roas(NW),rams(NW)) do 77 iw=1,NW rams(iw)=0.0d0 77 roas(iw)=0.0d0 allocate(Pt(ndim,ndim)) do 78 iw=1,ntt write(6,*)iw call tz(Pt,ndim,ndim) call MDE(n,ndim,x,Itr,Pt,ntw,tolw,ite,wl(iw)) call extr(Pt,alphar,gpr,ar,ndim,r,n) 78 call DORAMAS(NW,GAMMA,alphar,Gpr,Ar,5320.0d0,WMIN,DW, 1tempcm,lgauss,wl(iw),rams,roas) OPEN(41,FILE='rams.prn') OPEN(42,FILE='roas.prn') do 80 i=1,NW WRITE(41,3000)WMIN+DW*dble(i-1),rams(i) 80 WRITE(42,3000)WMIN+DW*dble(i-1),roas(i) 3000 FORMAT(f9.2,g12.4) close(41) close(42) end c ============================================================ SUBROUTINE readpol(a,im,fn,n) implicit none integer*4 im,n,i,j character*(*)fn real*8 a(n,3,3) open(90,file=fn) read(90,*) read(90,*) read(90,*)((a(im,i,j),i=1,j),j=1,3) close(90) do 1 j=1,3 do 1 i=1,j 1 a(im,j,i)=a(im,i,j) return end c ============================================================ SUBROUTINE readpolg(a,im,fn,n) implicit none integer*4 im,n,i,j character*(*)fn real*8 a(n,3,3) open(90,file=fn) read(90,*) read(90,*) read(90,*)((a(im,j,i),i=1,3),j=1,3) close(90) return end c ============================================================ SUBROUTINE readpola(a,im,fn,n) implicit none integer*4 im,n,i,j,k character*(*)fn real*8 a(n,3,3,3) open(90,file=fn) read(90,*) read(90,*) do 1 k=1,3 read(90,*)((a(im,k,i,j),i=1,j),j=1,3) do 1 j=1,3 do 1 i=1,j 1 a(im,k,j,i)=a(im,k,i,j) c first index electric close(90) return end c ============================================================ c ============================================================ SUBROUTINE READRX(r,n,fn,im) IMPLICIT none integer*4 Nat,n,i,im,j real*8 r(n,3),q,x,y,z,bohr CHARACTER*(*) fn(n) bohr=0.529177d0 OPEN(22,FILE=fn(im)) READ(22,*) read(22,*)Nat q=0.0d0 r(im,1)=0.0d0 r(im,2)=0.0d0 r(im,3)=0.0d0 do 1 i=1,Nat READ(22,*)j,x,y,z x=x/bohr y=y/bohr z=z/bohr q=q+dble(j) r(im,1)=r(im,1)+x*dble(j) r(im,2)=r(im,2)+y*dble(j) 1 r(im,3)=r(im,3)+z*dble(j) r(im,1)=r(im,1)/q r(im,2)=r(im,2)/q r(im,3)=r(im,3)/q write(6,600)im,Nat,(r(im,i),i=1,3) 600 format(' center ',i4,', ',i4,' atoms',3F9.4) close(22) RETURN END c ============================================================ subroutine report(s) character*(*) s write(6,*)s stop end c ============================================================ subroutine mktttt(rt,ttt) implicit none real*8 rt(3),ttt(3,3,3,3),r2,r4,r5,r9,delta integer*4 ix,jx,kx,lx r2=rt(1)**2+rt(2)**2+rt(3)**2 r4=r2*r2 r5=r4*dsqrt(r2) r9=r5*r4 do 613 ix=1,3 do 613 jx=1,3 do 613 kx=1,3 do 613 lx=1,3 613 ttt(ix,jx,kx,lx)=3.0d0*(35.0d0*rt(ix)*rt(jx)*rt(kx)*rt(lx) 1-5.0d0*r2*(delta(ix,kx)*rt(jx)*rt(lx) 1 +delta(ix,lx)*rt(jx)*rt(kx) 1 +delta(ix,jx)*rt(kx)*rt(lx) 1 +delta(kx,jx)*rt(lx)*rt(ix) 1 +delta(kx,lx)*rt(ix)*rt(jx) 1 +delta(jx,lx)*rt(ix)*rt(kx)) 1+r4*(delta(ix,lx)*delta(kx,jx) 1 +delta(kx,lx)*delta(ix,jx) 1 +delta(jx,lx)*delta(ix,kx)))/r9 return end c ============================================================ function delta(i,j) implicit none real*8 delta integer*4 i,j if(i.eq.j)then delta=1.0d0 else delta=0.0d0 endif return end c ============================================================ subroutine mkttt(rt,ttt) implicit none real*8 rt(3),ttt(3,3,3),rt2,rt5 integer*4 ix,jx,kx rt2=rt(1)**2+rt(2)**2+rt(3)**2 rt5=rt2**2*dsqrt(rt2) do 613 ix=1,3 do 613 jx=1,3 do 613 kx=1,3 ttt(ix,jx,kx)=-15.0d0*rt(ix)*rt(jx)*rt(kx)/rt5/rt2 if(ix.eq.jx)ttt(ix,jx,kx)=ttt(ix,jx,kx)+3.0d0*rt(kx)/rt5 if(ix.eq.kx)ttt(ix,jx,kx)=ttt(ix,jx,kx)+3.0d0*rt(jx)/rt5 613 if(jx.eq.kx)ttt(ix,jx,kx)=ttt(ix,jx,kx)+3.0d0*rt(ix)/rt5 return end c ============================================================ subroutine mktt(rt,tt) implicit none real*8 rt(3),tt(3,3),rt2,rt5 integer*4 ix,jx rt2=rt(1)**2+rt(2)**2+rt(3)**2 rt5=rt2**2*dsqrt(rt2) do 611 ix=1,3 do 611 jx=1,3 611 tt(ix,jx)=3.0d0*rt(ix)*rt(jx)/rt5 do 610 ix=1,3 610 tt(ix,ix)=tt(ix,ix)-rt2/rt5 return end c ============================================================ 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 DORAMAS(NW,GAMMA,alpha,G,A,EXCA,WMIN,DW,tempcm, 1lgauss,ECM,rams,roas) implicit none integer*4 NW,I real*8 alpha(3,3),G(3,3),A(3,3,3),EXCA,ECM,sf, 3roa1,ram1,WMIN,DW,tempcm,const,GAMMA,t real*8 rams(*),roas(*) logical lgauss call rr(ram1,roa1,alpha,G,A,EXCA) if(ECM.gt.1.0d0)then const=1.0d0/ECM/(1.0d0-exp(-ECM/tempcm)) do 4 i=1,NW t=sf(GAMMA,lgauss,WMIN+DW*dble(i-1),ECM)*const rams(i)=rams(i)+t*ram1 4 roas(i)=roas(i)+t*roa1 endif RETURN END c ============================================================ subroutine doX(n,r,X) implicit none integer*4 i,n,istart,j,jstart,ix,jx,kx,lx,ii,jj real*8 rij(3),r(n,3),tt(3,3),ttt(3,3,3),tttt(3,3,3,3), 1X(n*24,n*24) do 6 i=1,n istart=24*(i-1)+1 do 6 j=1,n jstart=24*(j-1)+1 if(i.ne.j)then do 61 ix=1,3 61 rij(ix)=r(i,ix)-r(j,ix) c T 0 0 0 -U 0 T=grad(1/rij), U=grad(T), G=grad(U), W=T/c c 0 T 0 0 0 -U c 0 0 W 0 0 0 c 0 0 0 W 0 0 c U 0 0 0 -G 0 c 0 U 0 0 0 -G call mktt(rij,tt) do 62 ix=1,3 do 62 jx=1,3 X(istart+ix-1, jstart+jx-1 )=tt(ix,jx) X(istart+ix-1+3,jstart+jx-1+3)=tt(ix,jx) X(istart+ix-1+6,jstart+jx-1+6)=tt(ix,jx)/137.0d0 62 X(istart+ix-1+9,jstart+jx-1+9)=tt(ix,jx)/137.0d0 call mkttt(rij,ttt) do 65 ix=1,3 ii=0 do 65 jx=1,3 do 65 kx=jx,3 ii=ii+1 X(istart+ix-1,jstart+11+ii)=-ttt(ix,jx,kx) X(istart+11+ii,jstart+ix-1)= ttt(ix,jx,kx) X(istart+ix+2,jstart+17+ii)=-ttt(ix,jx,kx) 65 X(istart+17+ii,jstart+ix+2)= ttt(ix,jx,kx) call mktttt(rij,tttt) ii=0 do 64 ix=1,3 do 64 jx=ix,3 ii=ii+1 jj=0 do 64 kx=1,3 do 64 lx=kx,3 jj=jj+1 X(istart+ii+11,jstart+11+jj)=-tttt(ix,jx,kx,lx) 64 X(istart+ii+17,jstart+17+jj)=-tttt(ix,jx,kx,lx) endif 6 continue return end c ============================================================ function sf(d,g,x,x0) implicit none logical g real*8 pi,spi,sf,d,x,x0,dd pi=4.0d0*atan(1.0d0) spi=dsqrt(pi) dd=((x-x0)/d)**2 c write(6,*)d,x,x0,g if(g)then if(dd.lt.20.0d0)then sf=exp(-dd)/d/spi else sf=0.0d0 endif else if(dd.lt.1.0d10)then sf=1.0d0/d/(dd+1.0d0)/pi else sf=0.0d0 endif endif return end c ============================================================ subroutine tz3(t,o,n,m) implicit none integer*4 o,n,m,i,j,k real*8 t(o,n,m) do 1 k=1,m do 1 j=1,n do 1 i=1,o 1 t(i,j,k)=0.0d0 return end c ============================================================ subroutine tz(t,n,m) implicit none integer*4 n,m,i,j real*8 t(n,m) do 1 j=1,m do 1 i=1,n 1 t(i,j)=0.0d0 return end c ============================================================ subroutine extr(Pt,alphar,gpr,ar,ndim,r,n) implicit none integer*4 i,j,k,l,ndim,n,kl,m,ii,pp,ip real*8 Pt(ndim,ndim),alphar(3,3),gpr(3,3), 2ar(3,3,3),rr(3),r(n,3),eps,delta,sum c call tz(alphar,3,3) call tz(gpr,3,3) call tz3(ar,3,3,3) c sum over all multipoles on particles i: do 8 i=1,n ii=24*(i-1) rr(1)=r(i,1) rr(2)=r(i,2) rr(3)=r(i,3) do 8 j=1,3 do 8 ip=1,n pp=24*(ip-1) kl=0 do 8 k=1,3 c from electric moment j = 1..3: c from electric field k = 1..3 on all particles ip: alphar(j,k)=alphar(j,k)+Pt(ii+j,pp+k) c from magnetic moment 7-9: c from electric field derivative 4..6: sum=0.0d0 do 85 l=1,3 do 85 m=1,3 85 sum=sum+eps(k,l,m)*rr(l)*Pt(ii+3+m,pp+3+j)/2.0d0 gpr(j,k)=gpr(j,k) -Pt(ii+6+k,pp+3+j)-sum do 8 l=k,3 kl=kl+1 c from quadrupole 13-18 c from electric field 1..3: c a 0 0 G A/3 0 c 0 a -G 0 0 A/3 c 0 -Gt 0 0 0 0 c Gt 0 0 0 0 0 c A/3 0 0 0 0 0 c 0 A/3 0 0 0 0 ar(j,k,l)=ar(j,k,l) 1 +3.0d0*Pt(ii+kl+12,pp+j) 1+1.5d0 *(rr(l)*Pt(ii+k ,pp+j) 1 +rr(k)*Pt(ii+l ,pp+j)) 1-delta(k,l)*(rr(1)*Pt(ii+1 ,pp+j) 1+ rr(2)*Pt(ii+2 ,pp+j) 1+ rr(3)*Pt(ii+3 ,pp+j)) 8 ar(j,l,k)=ar(j,k,l) return end c ============================================================ subroutine rr(ram1,roa1,alpha,G,A,EXCA) implicit none integer*4 I,J,K,L real*8 SAL0,SAL1,SAG0,SAG1,SA1,ram1,roa1,alpha(3,3), 1G(3,3),A(3,3,3),S180,D180,D,S90x,D90X,DX,S90Z,D90Z,DZ, 2YDY,YDX,AMU,BOHR,SpA3,beta2,gpisvejc,EXCA,betaa,betag,EPS SAL0=0.0d0 SAL1=0.0d0 SAG0=0.0d0 SAG1=0.0d0 SA1=0.0d0 DO 3 I=1,3 DO 3 J=1,3 SAL0=SAL0+alpha(I,I)*alpha(J,J) SAL1=SAL1+alpha(I,J)*alpha(I,J) SAG0=SAG0+alpha(I,I)* G(J,J) SAG1=SAG1+alpha(I,J)* G(I,J) DO 3 K=1,3 DO 3 L=1,3 3 SA1=SA1+EPS(I,K,L)*alpha(I,J)*A(K,L,J)/3.0d0 c c IL+IR, backscattering, DCP_I c SDCP180=3.0d0*SAL1-SAL0 c c IL+/-IR, backscattering S180=7.0d0*SAL1+SAL0 D180=8.0d0*(3.0d0*SAG1-SAG0+SA1) D=D180 IF(S180.gt.1.0d-9)D=D/S180 c c IL+/-IR, 90o, x S90X=7.0d0*SAL1+SAL0 D90X=2.0d0*(7.0d0*SAG1+SAG0+SA1) DX=D90X IF(S90X.gt.1.0d-9)DX=DX/S90X c c IL+/-IR, 90o, z S90Z=2.0d0*(3.0d0*SAL1-SAL0) D90Z=4.0d0*(3.0d0*SAG1-SAG0-SA1) DZ=D90Z IF(S90Z.gt.1.0d-9)DZ=DZ/S90Z c c polarized backscattering components (YDY+YDX=S180): YDY=3.0d0*SAL1-SAL0 YDX=4.0d0*SAL1+2.0d0*SAL0 AMU=1822.0d0 BOHR=0.529177d0 SpA3=dsqrt(SAL0)/3.0d0*dsqrt((AMU*BOHR**4)) beta2=0.5d0*(3*SAL1-SAL0)*(AMU*BOHR**4) YDY=6.0d0*YDY*(AMU*BOHR**4) YDX=6.0d0*YDX*(AMU*BOHR**4) c gpisvejc=980.0 gpisvejc=(AMU*BOHR**5)*1.0d4*2.0d0*4.0d0*atan(1.0d0)/EXCA c c alphag=SAG0/9.0d0*gpisvejc c c beta(A)^2=Delta2 in Gaussian: betaa=SA1*3.0d0/2.0d0*gpisvejc betag=(3.0d0*SAG1-SAG0)/2.0d0*gpisvejc c c DCP 180 c ram3=24.0d0*beta2 c roa3=8.0d0*(12.0d0*betag+4.0d0*betaa) c c unknown: S90X=S90X*(AMU*BOHR**4) D90X=D90x*gpisvejc DX=D90X IF(S90X.gt.1.0d-9)DX=DX/S90X c c ICP(90)=90Z c roa2=8.0d0*(3.0d0*betag-betaa) c ram2=12.0d0*beta2 c c ICP(180): c 8/c*(12*beta(G')^2 + 4*beta(A)^2): roa1=96.0d0*(betag+betaa/3.0d0) ram1=4.0d0*(45.0d0*SpA3**2+7.0d0*beta2) return end C ==========================================