SUBROUTINE ALPHA(AKR,CALPHA,CXALDS,CXALPH,CXALOF,CXALOS, & CXE0R,CXEPS,CXSC,CXSCR1,CXZC,CXZW,DX,IBETH,IBETH1, & ICOMP,IOCC,IPHI,IPHI1,IXYZ0,JORTH,MYID,MXCOMP,MXNAT,MXN3, & NAT,NAT0,NAT3,NCOMP,NX,NY,NZ,CXRLOC,CSHAPE,SHPAR) IMPLICIT NONE C*** Arguments: CHARACTER*6 CALPHA,CSHAPE INTEGER IBETH,IBETH1,IPHI,IPHI1,JORTH,MXCOMP,MXNAT,MXN3,MXMEM, & MYID,NAT,NAT0,NAT3,NCOMP,NX,NY,NZ COMPLEX CXALOF(NAT,3),CXALPH(NAT,3),CXE0R(3),CXEPS(MXCOMP), & CXALDS(NAT,3),CXALOS(NAT,3),CXSC(NAT,3,3),CXSCR1(NAT,3), & CXZC(NX+1,NY+1,NZ+1,6),CXZW(2*NX,2*NY,2*NZ,3), & CXRLOC(MXCOMP+1,3,3) INTEGER*2 ICOMP(NAT,3),IOCC(NAT),IXYZ0(MXNAT,3) REAL AKR(3),DX(3),SHPAR(6) C C*** Local Variables: LOGICAL INIT CHARACTER CMSGNM*70 INTEGER IA,IC,L,M,N,P,INFO,MINVIWRK(3) REAL A1(3),A2(3),A3(3),AK1,AK2,B1,B2,EMKD,PI,SUM, & SHAPEL(3,3),SHAPELT(3,3),WREV(3) COMPLEX CXI,CXRR,CXTERM,CX33(3,3), & CXALPH2(3,3),CXTMP(3,3),CXTMP2(3,3), & MINVZWRK(3) DOUBLE PRECISION ASHAP(3),LSHAP(3) C SAVE INIT,CXI,PI DATA CXI/(0.,1.)/ DATA INIT/.TRUE./ C*** Common: CHARACTER*6 CMETHD INTEGER DUM1,DUM2,DUM3,DUM4,DUM5,DUM6,DUM7,DUM8,DUM9, & DUM10,DUM11,DUM12 C----------------------------------------------------------------------- COMMON /M6/ DUM1, DUM2, DUM3, DUM4, DUM5, DUM6, DUM7, DUM8, & DUM9, DUM10, DUM11, DUM12, CMETHD C----------------------------------------------------------------------- C*********************************************************************** C Given: C AKR(1-3)=(kx,ky,kz)*d, where d=effective lattice spacing C CALPHA = 'LATTDR' or 'SCLDR' = polarizability prescription C CMETHD = determines 3-d FFT routine used by ESELF C CSHAPE = descriptor of target shape C CXE0R(1-3) = polarization vector in lattice coordinates C (assumed to be normalized) C CXEPS(1-NCOMP)=distinct values of dielectric constant C CXSC(1-NAT,1-3,1-3) = complex scratch space for SCLDR calculation C CXSCR1(1-NAT,1-3) = complex scratch space for SCLDR calculation C CXZC = complex scratch space needed by ESELF C CXZW = complex scratch space needed by ESELF C DX(1-3)=(dx/d,dy/d,dz/d), where dx,dy,dz=lattice spacings in C x,y,z directions, and d=(dx*dy*dz)**(1/3) C IBETH= MYID+IBETH1 if first time through combined BETA/THETA C orientation loop C IBETH1=starting value of IBETH (see above) C ICOMP(1-NAT3)=composition identifier for each lattice site and C direction (ICOMP=0 if lattice site is unoccupied) C storage scheme: C 1x,2x,...,NATx,1y,2y,...,NATy,1z,2z,...,NATz C IOCC(1-NAT) == 0 if site unoccupied C == 1 if site occupied C IPHI = IPHI1 if first time through PHI orientation loop C IPHI1= starting value of IPHI (see above) C JORTH= 1 if first incident polarization state C MXCOMP = dimensioning information C MXN3 = dimensioning information C MYID = parallel process identifier (=0 if only 1 process) C NAT = number of sites in extended target (incl. vacuum sites) C NAT3 = 3*number of sites in extended target (incl. vacuum sites) C NCOMP= number of different dielectric tensor elements in target C NX = x-dimension of extended target C NY = y-dimension of extended target C NZ = z-dimension of extended target C SHPAR(1-6)=target shape parameters C C Returns: C CXALPH(J,1-3)=(alpha_11,alpha_22,alpha_33)/d^3 for dipole J=1-NAT C where alpha_ij=complex polarizability tensor. C CXALOF(J,1-3)=(alpha_23,alpha_31,alpha_12)/d^3 for dipole J=1-NAT C C*** C If CALPHA = LATTDR: C Compute dipole polarizability using "Lattice Dispersion Relation" C of Draine & Goodman (1993,ApJ,March 10). It is required that C polarizability be such that an infinite lattice of such dipoles C reproduce the continuum dispersion relation for radiation C propagating with direction and polarization of radiation incident C on the DDA target. This is the recommended option. C C*** C Note: CXALPH = polarizability/d^3 C In the event that ICOMP=0, then we set CXALPH=1. C*********************************************************************** C B.T.Draine, Princeton Univ. Obs. C History: C 90.09.13 (BTD): Corrected error in DO loop limit. C 90.11.01 (BTD): Special treatment for "vacuum" sites in order C to allow FFT treatment. C 90.11.02 (BTD): Set CXALPH=(1.,0.) at vacuum sites. C 91.04.30 (BTD): Modified to include both O[(kd)^2] correction term C from Goedecke & Obrien (1988) in addition to C radiative reaction correction. C 91.05.07 (BTD): Modified to allow easy choice among three methods C for computing polarizability. C 91.05.07 (BTD): Added call to WRIMSG to record which method in use C 91.05.08 (BTD): Added CALPHA to argument list. C 91.05.09 (BTD): Added printing of kd and magnitude of correction C terms. C 91.05.13 (BTD): Experiment with use of approximate directional C average for coefficient of CXEPS in Lattice C Disperision Relation theory C 91.05.13 (BTD): Corrected error in radiative reaction correction C (error presumably introduced during last week) C 91.05.14 (BTD): Added separate options LDRXYZ and LDRAVG C 91.05.15 (BTD): Added option LDR000 to omit CXEPS-dependent C correction term C 91.05.30 (BTD): Changed LDRXYZ -> LDR100 C Added option LDR111 C 91.09.12 (BTD): Corrected error in numerical coefficient for C LDRAVG case C 91.09.17 (BTD): Eliminate options LDR000,LDRXYZ,LDRAVG C Introduce option LATTDR (LATTice Dispersion Relation) C which includes explicit dependence C on direction and polarization state C remove AK1 from argument list C add AKR, CXE0R to argument list C 92.05.14 (BTD): Introduce option LDRISO to experiment with C using average directional correction term C rather than using correction for incident direction C and polarization C 92.05.16 (BTD): Correct error in LDRISO option: average SUM should be C 3/15 rather than 4/15 C 97.11.02 (BTD): Add DX to argument list to allow use with noncubic C lattices (additional modifications still required!!) C 97.12.26 (BTD): Add CXALOF to argument list to allow use with noncubic C lattices (additional modifications still required). C 97.12.30 (BTD): Add code to evaluate alpha for noncubic case, C using subroutine NONCUBIC to evaluate the lattice C sums R0,R1,R2,R3 C 98.01.14 (BTD): Change sign of radiative reaction correction term. C 98.01.19 (BTD): Change treatment of R_1,C3, and C4 C and output value of C4 C temporary modification to allow value of C4 to be C read in from file 'alpha.par' C 98.01.22 (BTD): Minor change in computation of radiative reaction C correction. C 98.03.10 (BTD): set C4=-C3, disable reading C4 from 'alpha.par' C 98.04.27 (BTD): declared IC1,IC2,IC3 as integers. C 99.02.09 (BTD): experiment with change in expression used to calculate C off-diagonal elements of polarizability tensor C experiment with C4=-C3-R1 to "cancel" the contribution C of R1 to the off-diagonal terms C 03.01.29 (BTD): remove code to calculate alpha for noncubic lattice C (hold for release in future version) C end history C Copyright (C) 1993,1996,1997,1998,1999 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** PI=4.*ATAN(1.) AK2=AKR(1)*AKR(1)+AKR(2)*AKR(2)+AKR(3)*AKR(3) AK1=SQRT(AK2) CXRR=-(AK1*AK2/1.5)*CXI C*** EMKD = |m|*k_0*d , where |m|=refractive index C k_0 = wave vector in vacuo C d = lattice spacing EMKD=SQRT(CABS(CXEPS(1)))*AK1 C*********************************************************************** C*** Lattice dispersion relation C IF(CALPHA.EQ.'LATTDR')THEN WRITE(CMSGNM,FMT='(A,F6.4)') & 'Lattice dispersion relation for |m|k_0d=',EMKD CALL WRIMSG('ALPHA ',CMSGNM) C*** Compute sum (a_j*e_j)^2 , where a_j=unit propagation vector C e_j=unit polarization vector SUM=0.0 DO L=1,3 SUM=SUM+(AKR(L)*CABS(CXE0R(L)))**2 ENDDO SUM=SUM/AK2 B1=-1.8915316*AK2 B2=(0.1648469-1.7700004*SUM)*AK2 DO L=1,3 DO IA=1,NAT IC=ICOMP(IA,L) IF(IC.GT.0)THEN C*** First compute Clausius-Mossotti polarizability: CXTERM=(.75/PI)*(CXEPS(IC)-1.)/(CXEPS(IC)+2.) C*** Determine polarizability by requiring that infinite lattice of C dipoles have dipersion relation of continuum. CXTERM=CXTERM/(1.+CXTERM*(B1+CXEPS(IC)*B2)) C*** Radiative-reaction correction: CXALPH(IA,L)=CXTERM/(1.+CXTERM*CXRR) c set off-diagonal terms to zero CXALOF(IA,L)=0. ELSEIF(IC.EQ.0)THEN C To avoid divisions by zero, etc., set CXALPH=1 for vacuum sites. CXALPH(IA,L)=1. CXALOF(IA,L)=0. ENDIF ENDDO ENDDO RETURN ELSE WRITE(0,*)'Error: invalid option for subroutine ALPHA' STOP ENDIF END SUBROUTINE PETR(X,B,WRK,LDA,IPAR,SPAR,MATVEC,CMATVEC) IMPLICIT NONE C solves Ax=b C C Input: C B(LOCLEN) - RHS (doesn't change on output) C WRK(LDA, LOCLEN) - scratch array C IPAR - defines integer paramteres (see documentation of PIM) C SPAR - define real parameters C LDA - leading dimension of arrays C MATVEC -- name of procedure for computing y=Ax C CMATVEC -- name of procedure for computing y= conj(A') x C C Output: C X(LOCLEN) C C Reference: C Petravic,M and Kuo-Petravic,G. 1979, JCP, 32,263 C History: C Late 80' (Coded by Draine) C Changes to include "cprod" (PJF+BTD) C 95.06.01 (PJF) re-written to conform to PIM C 95.08.11 (BTD) minor editing (comments, etc.) C C Parameters: INTEGER IACE,IAXI,IGI,IPI,IQI,IR PARAMETER (IACE=1,IGI=2,IPI=3,IQI=4,IAXI=5,IR=6) C Arguments: INTEGER LDA INTEGER IPAR(*) REAL SPAR(*) COMPLEX B(*), WRK(LDA,*), X(*) EXTERNAL CMATVEC,MATVEC C C Local variables COMPLEX CDUMMY INTEGER IDUMMY,ITNO,L,LOCLEN,MAXIT,STATUS REAL & ALPHAI,BETAI1,BNORM,DUMMY,EPSILON,EXITNORM, & GIGI,GI1GI1,QIQI,RHSSTOP REAL SCSETRHSSTOP EXTERNAL PSCNRM2 C C*********************************************************************** LOCLEN = IPAR(4) MAXIT = IPAR(10) EPSILON=SPAR(1) C C Compute RHSSTOP = EPSILON*|B| C RHSSTOP = SCSETRHSSTOP(B,CDUMMY,EPSILON,IPAR,CDUMMY,PSCNRM2) C C Compute conjg(A')*B C CALL CMATVEC(B,WRK(1,IACE),IDUMMY) BNORM = 0. DO 170 L = 1,LOCLEN BNORM = BNORM + REAL(B(L)*CONJG(B(L))) WRK(L,IGI) = WRK(L,IACE) WRK(L,IPI) = WRK(L,IGI) 170 CONTINUE C C Compute |QI>=A|PI> C CALL MATVEC(WRK(1,IPI), WRK(1,IQI), IDUMMY) QIQI = 0. GIGI = 0. DO 190 L = 1,LOCLEN GIGI = GIGI + REAL(CONJG(WRK(L,IGI))*WRK(l,IGI)) QIQI = QIQI + REAL(CONJG(WRK(L,IQI))*WRK(L,IQI)) 190 CONTINUE ALPHAI = GIGI/QIQI DO 200 L = 1,LOCLEN X(L) = X(L) + ALPHAI*WRK(L,IPI) 200 CONTINUE C C compute |AX1>: C CALL MATVEC(X,WRK(1,IAXI),IDUMMY) C DO 340 ITNO = 2,MAXIT IPAR(11) = ITNO C C Transfer -> : C GI1GI1 = GIGI C C compute |GI>=AC|E>-AC|AXI>: C CALL CMATVEC(WRK(1,IAXI), WRK(1,IGI), IDUMMY) C C Compute GIGI=: GIGI = 0. DO 230 L = 1,LOCLEN WRK(L,IGI) = WRK(L,IACE) - WRK(L,IGI) 230 CONTINUE DO 235 L = 1,LOCLEN GIGI = GIGI + REAL(CONJG(WRK(L,IGI))*WRK(L,IGI)) 235 CONTINUE C C Compute BETAI-1=/: C BETAI1 = GIGI/GI1GI1 C C Compute |PI>=|GI>+BETAI-1|PI-1>: C DO 240 L = 1,LOCLEN WRK(L,IPI) = WRK(L,IGI) + BETAI1*WRK(L,IPI) 240 CONTINUE C C Compute |QI>=A|PI>: C CALL MATVEC(WRK(1,IPI), WRK(1,IQI), IDUMMY) C C Compute : C QIQI = 0. DO 250 L = 1,LOCLEN QIQI = QIQI + REAL(CONJG(WRK(L,IQI))*WRK(L,IQI)) 250 CONTINUE C C Compute ALPHAI = /: C ALPHAI = GIGI/QIQI C C Compute |XI+1>=|XI>+ALPHAI*|PI>: C DO 260 L = 1,LOCLEN X(L) = X(L) + ALPHAI*WRK(L,IPI) 260 CONTINUE C C Except every 10TH iteration, compute |AXI+1>=|AXI>+ALPHAI*|QI> C Draine: warning this part perhaps not checked C IF (ITNO.NE.10* (ITNO/10)) THEN DO 270 L = 1,LOCLEN WRK(L,IAXI) = WRK(L,IAXI) + ALPHAI*WRK(L,IQI) 270 CONTINUE ELSE CALL MATVEC(X,WRK(1,IAXI),IDUMMY) ENDIF C C Compute residual vector |RI>=A|XI>-|B> C DO 300 L=1, LOCLEN WRK(L,IR) = WRK(L,IAXI)-B(L) 300 CONTINUE C C Call STOPCRIT to check whether to terminate iteration. C Criterion: terminate if < TOL**2 * C Return with STATUS=0 if this condition is satisfied. C CALL STOPCRIT(B,WRK(1,IR),CDUMMY,CDUMMY,CDUMMY,WRK, + RHSSTOP,IDUMMY,EXITNORM, + STATUS,IPAR,DUMMY,DUMMY,DUMMY,DUMMY,PSCNRM2) C C Call PROGRESS to report "error" = sqrt()/sqrt() C CALL PROGRESS(LOCLEN,ITNO,EXITNORM,CDUMMY,CDUMMY,CDUMMY) C IF(STATUS.EQ.0)GOTO 500 340 CONTINUE 500 CONTINUE IPAR(12) = 0 RETURN END SUBROUTINE COPYIT(CXA,CXB,NPTS) IMPLICIT NONE C*** Arguments: INTEGER NPTS COMPLEX CXA(*), CXB(*) C C*** Local variables: INTEGER IPTS C*********************************************************************** C Purpose: To copy a complex vector from one memory location to another. C*********************************************************************** DO 10 IPTS=1,NPTS CXB(IPTS)=CXA(IPTS) 10 CONTINUE RETURN END SUBROUTINE DIAGL(U,V,IPAR) IMPLICIT NONE C Arguments: COMPLEX U(*), V(*) INTEGER IPAR(*) C Common variables: CHARACTER*6 CMETHD INTEGER & IDVOUT,MXN3F,MXNATF,MXNXF,MXNYF,MXNZF, & NAT,NAT0,NAT3,NX,NY,NZ COMPLEX CXADIA(1) C----------------------------------------------------------------------- COMMON/M2/CXADIA COMMON/M6/MXNATF,MXNXF,MXNYF,MXNZF,NAT,NAT3,NAT0,NX,NY,NZ, & MXN3F,IDVOUT,CMETHD C----------------------------------------------------------------------- C cg package C Left preconditioning subroutine - division by diagonal C History: C Originally written by P.J. Flatau, 1993 C 00.06.26 (BTD) Removed MXMEMF from COMMON/M6/ C end history c----------------------------------------------------------------------- CALL CDIV(U,V,CXADIA,NAT3) RETURN END SUBROUTINE CDIV(U,V,CXA,N) IMPLICIT NONE C Arguments: INTEGER N COMPLEX U(*), V(*), CXA(*) C C Local variables: INTEGER I c----------------------------------------------------------------------- c Part of conjugate gradient package c needed by left preconditioning subroutine c Copyright (C) 2003 B.T. Draine and P.J. Flatau c This code is covered by the GNU General Public License c----------------------------------------------------------------------- DO I=1,N V(I)=U(I)/CXA(I) ENDDO RETURN END SUBROUTINE MATVEC(CXX,CXY,IPAR) IMPLICIT NONE C Arguments: COMPLEX CXX(*), CXY(*), IPAR(*) C Common: INTEGER*2 IOCC(1) INTEGER & IDVOUT,MXN3F,MXNATF,MXNXF,MXNYF,MXNZF, & NAT,NAT0,NAT3,NX,NY,NZ CHARACTER*6 CMETHD REAL AK(3),DX(3) COMPLEX CXADIA(1),CXAOFF(1),CXZC(1),CXZW(1) C----------------------------------------------------------------------- COMMON/M1/AK,DX COMMON/M2/CXADIA COMMON/M3/CXZC COMMON/M4/CXZW COMMON/M5/IOCC COMMON/M6/MXNATF,MXNXF,MXNYF,MXNZF,NAT,NAT3,NAT0, & NX,NY,NZ,MXN3F,IDVOUT,CMETHD COMMON/M7/CXAOFF C----------------------------------------------------------------------- C Subroutine MATVEC C matvec --> cxy = A cxx 'N' C tmatvec --> cxy = A' cxx 'T' C cmatvec --> cxy = conjg(A') cxx 'C' C C C Notice that DDSCAT has cases 'N' and 'C' implemented in CPROD C Notice, however, that in DDSCAT 'T'='N'. C Some CG methods use cases N, C (like Petravic). C Some CG methods use cases N, T (like PIM) C Some CG use all three products N, C, T (parts of CCGPAK) C C Information about matrix A is transfered via commons /M1/-/M7/ with C the main program. Arrays cxadia, cxzc, cxzwm, iocc are properly C dimensioned inside "cprod". You may get warning during the compilation C about "misaligned common". Ignore these. Also notice mxnatf, C mxn3f,mxnxf, mxnyf, mxnzf. These are defined by "parameter C statement" in main code and this is the reason why they have "f" suffix. C C ipar - information array, not used, defined for compatibility C C History: C 97.11.01 (BTD): modified to allow for noncubic lattice, introduced C new vector DX(1-3) via COMMON/M1/ C added DX to argument list of CPROD C 97.12.25 (BTD): added COMMON/M7/CXAOFF to communicate off-diagonal C elements of inverse polarizability tensors. C Added CXAOFF to argument list of CPROD C 00.06.13 (BTD): Explicit declaration of all variables. C 00.06.25 (BTD): Removed MXMEMF from COMMON/M6/ C end history C Copyright (C) 1993,1997,2000 C B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** CALL CPROD(AK,DX,CMETHD,'N',CXADIA,CXAOFF, & CXX,CXY,CXZC,CXZW, & IDVOUT,IOCC,MXN3f,MXNATf,MXNXf,MXNYf,MXNZf, & NAT,NAT0,NAT3,NX,NY,NZ) RETURN END SUBROUTINE CMATVEC(CXX, CXY, IPAR) IMPLICIT NONE C Arguments: COMPLEX CXX(*), CXY(*), IPAR(*) C Common: INTEGER*2 IOCC(1) INTEGER & IDVOUT,MXN3F,MXNATF,MXNXF,MXNYF,MXNZF, & NAT,NAT0,NAT3,NX,NY,NZ CHARACTER*6 CMETHD REAL AK(3),DX(3) COMPLEX CXADIA(1),CXAOFF(1),CXZC(1), CXZW(1) C----------------------------------------------------------------------- COMMON/M1/AK,DX COMMON/M2/CXADIA COMMON/M3/CXZC COMMON/M4/CXZW COMMON/M5/IOCC COMMON/M6/MXNATF,MXNXF,MXNYF,MXNZF,NAT,NAT3,NAT0,NX,NY,NZ, & MXN3F,IDVOUT,CMETHD COMMON/M7/CXAOFF C----------------------------------------------------------------------- c Subroutine CMATVEC c C matvec --> cxy = A cxx C tmatvec --> cxy = A' cxx C cmatvec --> cxy = conjg(A') cxx C History: C 97.11.01 (BTD): modified to allow for noncubic lattice, introduced C new vector DX(1-3) via COMMON/M1/ C added DX to argument list of CPROD C 97.12.25 (BTD): added COMMON/M7/CXAOFF to communicate off-diagonal C elements of inverse polarizability tensors. C Added CXAOFF to argument list of CPROD C 00.06.13 (BTD): Explicit declaration of all variables. C 00.06.26 (BTD): Removed MXMEMF from COMMON/M6/ C end history C Copyright (C) 1993,1997,2000 C B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** CALL CPROD(AK,DX,CMETHD,'C',CXADIA,CXAOFF, & CXX,CXY,CXZC,CXZW, & IDVOUT,IOCC,MXN3F,MXNATF,MXNXF,MXNYF,MXNZF, & NAT,NAT0,NAT3,NX,NY,NZ) RETURN END SUBROUTINE CPROD(AK,DX,CMETHD,CWHAT,CXADIA,CXAOFF, & CXX,CXY,CXZC,CXZW, & IDVOUT,IOCC,MXN3,MXNAT,MXNX,MXNY,MXNZ, & NAT,NAT0,NAT3,NX,NY,NZ) IMPLICIT NONE c Arguments: INTEGER IDVOUT,MXN3,MXNAT,MXNX,MXNY,MXNZ, & NAT,NAT0,NAT3,NX,NY,NZ CHARACTER CMETHD*6,CWHAT*1 INTEGER*2 IOCC(MXNAT) REAL AK(3),DX(3) COMPLEX CXADIA(MXN3),CXAOFF(MXN3), & CXX(MXN3),CXY(MXN3),CXZC(NX+1,NY+1,NZ+1,6), & CXZW(MXNX*MXNY*MXNZ,*) c Local variables: INTEGER J1,J2,J3 REAL AKD C Intrinsic functions: INTRINSIC CONJG C*********************************************************************** c Subroutine CPROD C Given: C AK(1-3)= k(1-3)*d , where k=k vector in vacuo C d=effective lattice spacing=(dx*dy*dz)**(1/3) C DX(1-3)=(dx/d, dy/d, dz/d) C where dx,dy,dz=lattice spacing in x,y,z directions. C By definition of effective lattice spacing d, must C have DX(1)*DX(2)*DX(3)=1 C CMETHD = character variable specifying method used for FFT evaluation C (see ESELF) C CWHAT = 'N' or 'C' determines what is to be done (see below) C CXADIA(1-NAT3) = diagonal elements of A(j,k) C CXAOFF(1-NAT3) = off-diagonal elements of 3x3 diagonal blocks of C A(j,k): CXAOFF(J,1-3)=a_{23},a_{31},a_{12} from 3x3 C matrix a_{ij} for dipole J. C Recall that a_{ij} is inverse of 3x3 polarizability C tensor alpha_{ij} for dipole J. C CXX(1-NAT3) = complex vector C CXZC = (NX+1)*(NY+1)*(NZ+1)*6 array of Green function C coefficients used internally by ESELF, and generated C by ESELF each time called with new k vector C *****Not to be overwritten between calls to ESELF***** C CXZW = complex workspace required by ESELF C IDVOUT = device number for output C IOCC(1-NAT) = 0 or 1 if site is unoccupied or occupied C MXN3,MXNX,MXNY,MXNZ = dimensioning information (see DDSCAT) C NAT = number of sites in extended target C NAT0 = number of occupied sites in extended target C NAT3 = 3*NAT C NX,NY,NZ = extended target size is (NX*DX(1)) by (NY*DX(2)) by C (NZ*DX(3)) C C If CWHAT = 'N', this returns CXY(j)=A(j,k)*CXX(k) (sum over k) C = 'C', this returns CXY(j)=conjg(A(k,j))*CXX(k) (sum over k) C C It is assumed that matrix A(j,k) is symmetric C Diagonal elements of A(j,k) are in vector CXADIA(j) C Off-diagonal terms coupling dipole component x with direction y, etc C are stored in vector CXAOFF(j), CXAOFF(j+NAT), CXAOFF(j+2*NAT) C Other off-diagonal elements of A(j,k) are produced internally by C subroutine ESELF, which uses 3d-FFT to accomplish convolution C of CXX with A-matrix terms coupling dipole j with other dipoles C (i.e., E field at j due to other dipoles) C C Original version of CPROD created by P.J.Flatau, Colorado State Univ. C Subsequently modified by B.T.Draine, Princeton Univ. Obs. C History: C 90.11.04 (BTD): Major modification to incorporate FFT method (FFT C computations are handled by ESELF). C 90.11.08 (BTD): Remove DO...ENDDO structures; renumber DO...CONTINUEs C 90.11.21 (BTD): Remove "CCGMVM" option: only FFT stuff retained. C 90.12.03 (BTD): No longer necessary to change ordering of CXX and CXY C 90.12.05 (BTD): Changed call to ESELF to eliminate MXNX,MXNY,MXNZ C 97.11.01 (BTD): Added comments and modified to accomodate noncubic C lattice. Added DX to argument list for CPROD and to C argument list in call to ESELF C 97.12.25 (BTD): Added CXAOFF to argument list. C 97.12.28 (BTD): Modified to use CXAOFF in calculation of CXY C 98.01.01 (BTD): Modified coding to use indices J1,J2,J3 in loops C computing with CXAOFF C 00.06.13 (BTD): Explicit declaration of all variables. C end history C Copyright (C) 1993,1997,1998,2000 C B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. c----------------------------------------------------------------------- AKD=SQRT(AK(1)*AK(1)+AK(2)*AK(2)+AK(3)*AK(3)) C We assume that input vector CXX has zeroes for elements corresponding C to vacuum sites. IF(CWHAT.EQ.'N')THEN C*********************************************************************** C CALL ESELF(CMETHD,CXX,NX,NY,NZ,AKD,DX,CXZC,CXZW,CXY) C C*********************************************************************** C ESELF omitted contribution from diagonal terms: add these. C Also remember that ESELF computes -A_jk*x_k DO J1=1,NAT3 CXY(J1)=CXADIA(J1)*CXX(J1)-CXY(J1) ENDDO C Add contribution from off-diagonal elements of 3x3 diagonal blocks C Version below assumes that off-diagonal elements a_ij are stored as C CXAOFF(J,1-3)=(a_23,a_31,a_12) for dipole J. C 97.12.28 (BTD): DO J1=1,NAT J2=J1+NAT J3=J2+NAT CXY(J1)=CXY(J1)+CXAOFF(J2)*CXX(J3)+CXAOFF(J3)*CXX(J2) CXY(J2)=CXY(J2)+CXAOFF(J3)*CXX(J1)+CXAOFF(J1)*CXX(J3) CXY(J3)=CXY(J3)+CXAOFF(J1)*CXX(J2)+CXAOFF(J2)*CXX(J1) ENDDO ELSEIF(CWHAT.EQ.'C')THEN C Need to temporarily replace CXX by conjg(CXX) DO J1=1,NAT3 CXX(J1)=CONJG(CXX(J1)) ENDDO C********************************************************************** C CALL ESELF(CMETHD,CXX,NX,NY,NZ,AKD,DX,CXZC,CXZW,CXY) C C************************************************************************ C C ESELF omitted contribution from diagonal terms: add these. C Also remember that ESELF computed -A_jk*conjg(x_k), while C we seek to compute conjg(A_kj)* x_k. Since A_kj=A_jk, we have C conjg(A_kj)*x_k=conjg(A_jk*conjg(x_k)) C DO J1=1,NAT3 CXY(J1)=CXADIA(J1)*CXX(J1)-CXY(J1) ENDDO C C Add contribution from off-diagonal elements of 3x3 diagonal blocks C It is assumed that we have stored off-diagonal elements a_ij as C CXAOFF(J,1-3)=(a_23,a_31,a_12) for element J C Note that we could *probably* speed up the code by C (1) "reducing" product vector CXY immediately after calling ESELF C (2) restricting CXADIA to occupied sites C (3) restricting CXAOFF to occupied sites C but this would require modifications outside this routine, and C therefore requires further study... C 97.12.28 (BTD) DO J1=1,NAT J2=J1+NAT J3=J2+NAT CXY(J1)=CXY(J1)+CXAOFF(J2)*CXX(J3)+CXAOFF(J3)*CXX(J2) CXY(J2)=CXY(J2)+CXAOFF(J3)*CXX(J1)+CXAOFF(J1)*CXX(J3) CXY(J3)=CXY(J3)+CXAOFF(J1)*CXX(J2)+CXAOFF(J2)*CXX(J1) ENDDO C Now compute conjugate of CXY, and restore CXX: DO J1=1,NAT3 CXY(J1)=CONJG(CXY(J1)) CXX(J1)=CONJG(CXX(J1)) ENDDO ELSE WRITE(IDVOUT,FMT=*)' Error in CPROD: CWHAT= ',CWHAT ENDIF C*** Now must "clean" product vector CXY: C (zero out elements corresponding to vacuum sites) IF(NAT0.LT.NAT)CALL NULLER(CXY,IOCC,MXNAT,MXN3,NAT) RETURN END SUBROUTINE CXC3DFFT(CZ,I1,I2,I3,I4,I5,I6,I7) IMPLICIT NONE COMPLEX CZ(1) INTEGER I1,I2,I3,I4,I5,I6,I7 C*********************************************************************** C Note: This module should NOT be linked into DDSCAT when C running on a Convex with the vector library -- in that case C link with the -lveclib flag and use the native Convex C3DFFT C routine. C On other systems, however, you need to link in this routine C to avoid an error message from the compiler. C C History: C 93.xx.xx (BTD) created to substitute for convex C3DFFT library routine C 00.07.05 (BTD) renamed for greater Makefile consistency C end history C Copyright (C) 1993,2000 C B.T. Draine and P.J. Flatau C C This code is covered by the GNU General Public License. C*********************************************************************** CALL ERRMSG('FATAL','C3DFFT',' CONVEX option but dummy routine!') RETURN END SUBROUTINE CXFFTW(CX,MX,MY,MZ,ISIGN) IMPLICIT NONE C C Arguments: C INTEGER ISIGN,MX,MY,MZ COMPLEX CX(MX,MY,MZ) C*********************************************************************** C C Purpose: This is a dummy routine to substitute for CXFFTW to allow C DDSCAT to be used on systems where the FFTW package has not C been installed. If called, it will generate a fatal error C with an error message reporting DDSCAT has been compiled C with this dummy routine and option FFTW is unavailable. C C History: C 00.07.05 (BTD): created C 03.07.13 (BTD): changed option FFTWFJ to option FFTW21 C end history C Copyright (C) 2000, 2003 C B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C C*********************************************************************** C WRITE(0,*)'FATAL ERROR: DDSCAT has been compiled with dummy', & ' version of CXFFTW' WRITE(0,*)' *** option FFTW21 cannot be used with dummy routine' WRITE(0,*)' *** to enable option FFTW21 it is necessary to:' WRITE(0,*)' --- have fftw 2.1.x library installed on system' WRITE(0,*)' --- in Makefile, variable LIBFFTW should include', & ' correct path to libsfftw.a' WRITE(0,*)' --- in Makefile, define cxfftw = cxfftw', & ' (rather than dummycxfftw.f)' WRITE(0,*)' --- make ddscat' STOP END SUBROUTINE WRITENET(CNETFILE,IDNC,CNETFLAG,CSTAMP,CDESCR,CMDFFT, & CALPHA,CSHAPE,CMDTRQ,MXWAV,MXRAD,MXTHET, & MXBETA,MXPHI,MXN3,MXNAT,MXSCA,MXTIMERS,IORTH, & NX,NY,NZ,NAT0,NAT,NAT3,NSCAT,NCOMP,NORI,NWAV, & NRAD,NTHETA,NPHI,NBETA,NTIMERS,MXNX,MXNY,MXNZ, & MXCOMP,BETMID,BETMXD,THTMID,THTMXD, & PHIMID,PHIMXD,ICTHM,IPHIM,TOL,DX,WAVEA,AEFFA, & THETA,BETA,PHI,A1,A2,ICOMP,IXYZ0,IWAV,IRAD, & ITHETA,IBETA,IPHI,IORI,AEFF,WAVE,XX,AK1,BETAD, & THETAD,PHID,AKR,CXE01R,CXE02R,CXRFR,CXEPS, & QEXT,QABS,QSCAT,G,QBKSCA,QPHA,QAV,QTRQAB, & QTRQSC,SM,SMORI,PHIN,THETAN,TIMERS) C*********************************************************************** C Dummy version of subroutine WRITENET C (genuine subroutine WRITENET writes NetCDF portable binary file C to file cnetfile 'dd.cdf') C C This dummy version is provided in order that use of DDSCAT is C possible even if NetCDF software is not installed on system. C C History: C 96.11.15 (PJF): First version of NetCDF (Wolff's request) C 5a6 ---> 5a7 Converted to subroutine. C 96.11.20 (BTD): dummy version created from real WRITENET C 98.01.11 (BTD): added DX(3) to argument list for compatibility with C new version of WRITENET C end history C Copyright (C) 1996,1998 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C C .. Scalar Arguments .. REAL AEFF,AK1,BETAD,BETMID,BETMXD,PHID,PHIMID,PHIMXD,THETAD, & THTMID,THTMXD,TOL,WAVE,XX INTEGER IBETA,ICTHM,IDNC,IORI,IORTH,IPHI,IPHIM,IRAD,ITHETA,IWAV, & MXBETA,MXCOMP,MXN3,MXNAT,MXNX,MXNY,MXNZ,MXPHI,MXRAD,MXSCA, & MXTHET,MXTIMERS,MXWAV,NAT,NAT0,NAT3,NBETA,NCOMP,NORI,NPHI, & NRAD,NSCAT,NTHETA,NTIMERS,NWAV,NX,NY,NZ CHARACTER CALPHA*6,CMDFFT*6,CMDTRQ*6,CNETFLAG*6,CSHAPE*6, & CNETFILE*14,CSTAMP*22,CDESCR*67 C .. Array Arguments .. COMPLEX CXE01R(3),CXE02R(3),CXEPS(MXCOMP),CXRFR(MXCOMP) REAL A1(3),A2(3),AEFFA(MXRAD),AKR(3),BETA(MXBETA),DX(3),G(2), & PHI(MXPHI),PHIN(MXSCA), & QABS(2),QAV(17),QBKSCA(2),QEXT(2),QPHA(2), & QSCAT(2),QTRQAB(3,2),QTRQSC(3,2),SM(MXSCA,4,4), & SMORI(MXSCA,4,4),THETA(MXTHET),THETAN(MXSCA), & TIMERS(MXTIMERS),WAVEA(MXWAV) INTEGER*2 ICOMP(MXN3),IXYZ0(MXNAT,3) C C This subroutine should never be called, so if it is, bail out... C CALL ERRMSG('FATAL','DUMMYWRITENET', & ' NetCDF is not installed. Change to NOTCDF in ddscat.par!') RETURN END SUBROUTINE TIMEIT(CMSGTM,DTIME) C C timeit_null C C This version of timeit is a dummy which does not provide any C timing information. C 94.06.20 (PJF) add dtime to the formal parameters C 99.05.12 (BTD) explicitly declare DTIME C C Copyright (C) 1993,1999 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C Arguments: CHARACTER CMSGTM*(*) REAL DTIME C Local variables: CHARACTER CSTA*3,CMSGNM*70 C*********************************************************************** DATA CSTA/'ONE'/ IF(CSTA.EQ.'ONE')THEN CSTA='TWO' ELSEIF(CSTA.EQ.'TWO')THEN CSTA='ONE' WRITE(CMSGNM,FMT='(A,A)')'Timing results for: ',CMSGTM CALL WRIMSG('TIMEIT',CMSGNM) IF(CMSGTM(1:7).NE.'NOPRINT'.AND.CMSGTM(1:7).NE.'noprint')THEN WRITE(CMSGNM,FMT='(A)') & 'Dummy timing routine: no timings available' CALL WRIMSG('TIMEIT',CMSGNM) ENDIF ENDIF RETURN END