gulp af
###links http://fortran90.org/src/best-practices.html#modules-and-programs
http://www.csee.umbc.edu/~squire/fortranclass/summary.shtml#Stru
#QUESTIONS: 01/06/12 ##Q
! use constants,
use current
use general, only : bfactor
use iochannels
use parallel
implicit none###question
where are the modules at?
####answer
Refer to file:
~/Src/modules.f90
##Q
subroutine thermalconductivity(mcv,derv2,eigr,Sij,freq,nphonatc,ncfoc,nphonatptr,maxd2)###question why is Sij listed with the subroutine definition, when it seems to me to be more of a "local" variable since its contents are defined within the thermalconductivtiy.f90? ####answer It could be local as things stand. However Sij is the result of the subroutine and so it is conceivable that it might be needed as a return argument. For example, if fitting of Sij was allowed then this might be the case. There is also the possibility that memory overlays can be used (i.e. an array that is used elsewhere in the calling routine "phonon" could be passed as the Sij workspace). ##Q
!
! Passed variables
!
integer(i4), intent(in) :: nphonatptr(*)
integer(i4), intent(in) :: maxd2
integer(i4), intent(in) :: mcv
integer(i4), intent(in) :: ncfoc
integer(i4), intent(in) :: nphonatc
real(dp), intent(inout) :: derv2(maxd2,*)
real(dp), intent(in) :: eigr(maxd2,*)
real(dp), intent(out) :: Sij(mcv,*)
real(dp), intent(in) :: freq(*)###question where do these variables get passed in from the main part of the code? this is one of the challenges for me, tracing functionality back to the part of the code which reads in the gulp.gin. ####answer These get passed in from "phonon" as the calling master phonon calculation routine. The best thing is rather than have you track things back in the code from the input file that I tell you what the variables are in a physical sense:
maxd2 - this is the left-hand dimension of the arrays derv2 and eigr
mcv - this is the left-hand dimension of Sij and is 3 x ncfoc
ncfoc - this is the number of fully occupied cores in the system. In the simplest case this is just equal to the number of atoms in the unit cell (numat). In cases where there is a shell model being used or there are partially occupied sites then it can be < numat.
nphonatc - this is the number of cores in the system. Again, in the simplest case it's just equal to numat.
nphonatptr(*) - this is an array of integers that point from the number of the core (from 1 to nphonatc) - to the number of full coordinate list site (that includes shells). For the simple case, then nphonatptr(i) = i.
derv2(maxd2,) - this is the dynamical matrix (i.e. mass weighted second derivatives) as a 3ncfoc x 3*ncfoc matrix
eigr(maxd2,) - these are the eigenvectors as a 3ncfoc x 3*ncfoc matrix
freq() - this is a linear array containing the frequencies in ascending order of the 3ncfoc modes
###question what is the meaning of the syntax "derv2(maxd2,*)"? what will the size and shape of the array be? ####answer
I think I answered my own question, line 164:
!
! Copy results back to Sij
!
Sij(1:mcv,1:mcv) = derv2(1:mcv,1:mcv)what is the variable "mcv"? It must be the number of modes in the system, 3num_atoms_total. Based on this, how do I know to label the passed array as "derv2(maxd2,)"? ##Q mcv is indeed 3 x num_atoms_total (numat in GULP). Not sure quite what the second question is asking?
!
! Allocate local array ldone to avoid duplicate multiplies in case of partial occupancy
!
ncfoc2 = ncfoc*(ncfoc+1)/2
allocate(ldone(ncfoc2),stat=status)
if (status/=0) call outofmemory('thermalconductivity','ldone')
ldone(1:ncfoc2) = .false.###question what does "partial occupancy" mean? I guess ldone(ncfoc2) is an array of logical ".true." or ".false." of length ncfoc2, based on whether the array was allocated correctly? ####answer Because 2 or more atoms in the full cell can be linked as being part of the same fully occupied site (due to partial occupancies) then their matrix elements have to be combined. The flag ldone indicates whether this site has been addressed already or not. ##Q
!
! Scale dynamical matrix elements by minimum image nearest distance between sites
!
do i = 1,nphonatc
ii = nphonatptr(i)
ix = 3*ii - 2
iy = ix + 1
iz = ix + 2
do j = 1,i-1
jj = nphonatptr(j)
ind = ii*(ii-1)/2 + jj
if (.not.ldone(ind)) then
jx = 3*jj - 2
jy = jx + 1
jz = jx + 2###question what is the meaning of the variable "nphonatc" and "nphonatptr"? I guess it they are related to either the number of modes, or number of atoms, or both. ####answer This is covered in an earlier answer. ##Q
!
! Compute initial vector
!
xd = xclat(jj) - xclat(ii)
yd = yclat(jj) - yclat(ii)
zd = zclat(jj) - zclat(ii)###question the variables "{x,y,z}clat" are defined as local variables, but don't seem to be allocated or initialized anywhere before these lines? Does this mean all values are = 0? ####answer {x,y,z}clat are actually global variables that come from modules. They are the x, y and z, respectively, Cartesian coordinates of the atoms in Angstroms.
#NEW
#QUESTIONS: 01/08/12
##Q
!
! Create inverse frequency factors while trapping translations and imaginary modes
!
allocate(freqinv(mcv),stat=status)
if (status/=0) call outofmemory('thermalconductivity','freqinv')
nfreqmin = 0
do i = 1,mcv
if (freq(i).gt.1.0_dp) then
if (nfreqmin.eq.0) nfreqmin = i
freqinv(i) = 1.0_dp/sqrt(2.0_dp*freq(i))
else
freqinv(i) = 1.0_dp
endif
enddo###question do you store imag freq as negative numbers also internally in the code? I know neg freq == imag in the output.gout. Also, why do you check "freq(i).gt.1.0_dp" and not "freq(i).gt.0.0_dp"? ####answer GULP does indeed store imaginary frequencies as a negative number for convenience. The reason for the check that the frequency > 1 is to eliminate the divide by zero error that would occur for the translational modes in the lines that follow. This might be a bit crude, but is unlikely that most systems would have a genuine frequency that was below 1 cm^-1. ##Q ###question What are the units of "freq" internally?
####answer The units of freq are in wave numbers (cm^-1) on return from pdiag(g). ##Q ###question could you add an option to output the eigenvectors like the freq:
output freq freq.gout
output eig eig.gout ####answer If you add the keyword "eigen" to the input file then you should get a print out of the eigenvectors in the .gout file. If you want them in a separate file then I could look to add this as an option, but would obviously take longer.
##Q ###question And also an option to output the D_i:
output Di Di.gout ####answer The Di values are given in the output of the program. However, I can again add an option to write to a separate file if you let me know what would be the preferred format.
##Q
###question
####answer
#result
##GULP 4.0.5
##GULP 4.0.7-units
###thermalconductivity_gale_022413
source code:
###thermalconductivity
I changed your defn of the lorentzian, which seems to be mising a factor of pi, to:
dwij = (1.0/pi)*(bfactor*dwavg)/( (freq(j) - freq(i))**2 + (bfactor*dwavg)**2 )and this result looks identical to what my matlab code gets, except for the unit conversion.
https://raw.github.com/jasonlarkin/disorder/master/matlab/gulp/lj_amor_4.0.7_units-lor.jpg
source code:
here is the input.gin for that structure if you want to play with it to figure out the units:
https://raw.github.com/jasonlarkin/disorder/master/matlab/gulp/amor_cond.gin
##correct
###units obviously the unit conversion above is no good, here are my factors:
function Di = m_gulp_af_lj_readDi(path,name)
con = m_constant; lj = m_lj;
Di = load([path name]);
%gets to m^2/s
factor.gulp = con.eV2J^2*con.avog^2/con.ang2m^2/con.c
factor =
gulp: 3.105287902711237e+19
con.eV2J
ans =
1.602176460000000e-19
con.avog
ans =
6.022141290000000e+23
con.ang2m
ans =
1.000000000000000e-10
con.c
ans =
2.997924580000019e+10
%gets to lj units
factor.lj = lj.tau/lj.sigma^2
Di_convert(:,3) = Di(:,3)*factor.gulp*factor.lj;
Di_convert(:,2) = Di(:,2)*con.c*lj.tau;
loglog(Di_convert(:,2),Di_convert(:,3),'.')
end##GULP 4.0.7-lor
I changed the following:
real(dp) :: pi
pi = 3.14159265359!
constant = 1.0d10*pi/12.0_dp ! 1/3 convoluted with 1/2 squared from A7!
! Scale by constants and inverse frequency squared - factor of third is for averaging over directions
!
Di(i) = Di(i) + Di_loc*constant/(freq(i)**2)This is alot closer:
#conv
The updated version of the thermalconductivity.f90:
https://github.com/jasonlarkin/disorder/blob/master/matlab/gulp/thermalconductivity.f90
I changed the broadening to the average freq spacing, and now the variable broaden controls what multiple of dwavg the Lor is broadened by:
! DEBUG
!
! Find mean level spacing
!
dwavg = 0.0_dp;
do i = nfreqmin,mcv-1
dwavg = freq(i+1) - freq(i) + dwavg
enddo
dwavg = dwavg/(mcv-1 - nfreqmin)
! DEBUG
!
! Sum over coupling with mode j weighted by Lorentzian factor
!
do j = nfreqmin,i-1
dwij = (bfactor*dwavg)/(1 + ((bfactor*dwavg)*(freq(j) - freq(i)))**2)
Di = Di + dwij*Sij(j,i)**2
enddo
do j = i+1,mcv
dwij = (bfactor*dwavg)/(1 + ((bfactor*dwavg)*(freq(j) - freq(i)))**2)
Di = Di + dwij*Sij(j,i)**2
enddo
I also added a print out of the values of i,j,rij:
! DEBUG
!
! Write out a file with rij
!
open(91,file='rij.gout',status='unknown',form='formatted')
open(91,file='rij.gout',status='unknown',form='formatted')
write(91,'(i6,2x,i6,2x,f12.4)') &
&i,j,rij
! DEBUGHere is a comparison for the following structure:
4 4 20.000 20.000 20.000
1 1 0.0000 0.0000 0.0000
2 1 0.7818 0.7818 0.0000
3 1 0.7818 0.0000 0.7818
4 1 0.0000 0.7818 0.7818##GULP
The results are in this folder:
https://github.com/jasonlarkin/disorder/tree/master/matlab/gulp/conv
##matlab
Here is the matlab code used to produce the below:
https://github.com/jasonlarkin/disorder/blob/master/matlab/m_af_lj.m
###freq
0+2.0197i
0+2.0197i
0+2.0197i
0+7.6649e-08i
1.0875e-07
2.9761e-07
8.6421
8.6421
12.387
12.387
12.387
17.635####gulp
-5.003599416318
-5.003599416318
-5.003599416318
-0.000000401238
-0.000000215338
0.000000055440
21.410294901138
21.410294901138
30.689370515317
30.689370515317
30.689370515317
43.688796417065which converts to lj units as:
-2.0195
-2.0195
-2.0195
-0.0000
-0.0000
0.0000
8.6414
8.6414
12.3866
12.3866
12.3866
17.6333###Di
0+2.0197i -0
0+2.0197i -0
0+2.0197i -0
0+7.6649e-08i -9.5417e-16-6.604e-30i
1.0875e-07 1.5531e-15-3.7303e-30i
2.9761e-07 1.1245e-16-3.4781e-31i
8.6421 0.052754
8.6421 0.052754
12.387 0.092533
12.387 0.092533
12.387 0.092533
17.635 0.11164###eigvec
For the eigenvectors, I suppose there is degeneracy (or near degeneracy), so the eigenvectors don't match up exactly:
####matlab
0.45418 0.20003 0.060894 0.47816 -0.0037812 -0.14611 -0.20412 0.35355 0.28183 0.35721 0.2073 -0.28868
-0.36588 0.26938 -0.20873 0.14424 -0.068707 0.47379 -0.20412 -0.35355 0.28336 0.2939 -0.28867 -0.28868
-0.088302 -0.46941 0.14783 -0.023661 -0.49524 -0.064615 0.40825 -3.4235e-16 -0.12093 0.4812 -0.061803 -0.28868
-0.20219 0.38301 0.24985 0.47816 -0.0037812 -0.14611 0.20412 -0.35355 -0.40429 0.18731 0.22687 0.28868
0.29049 0.086404 -0.39768 0.14424 -0.068707 0.47379 0.20412 0.35355 -0.40275 0.12399 -0.2691 0.28868
0.088302 0.46941 -0.14783 -0.023661 -0.49524 -0.064615 0.40825 -6.1542e-16 0.12093 -0.4812 0.061803 -0.28868
0.20219 -0.38301 -0.24985 0.47816 -0.0037812 -0.14611 0.20412 -0.35355 0.40429 -0.18731 -0.22687 0.28868
0.36588 -0.26938 0.20873 0.14424 -0.068707 0.47379 -0.20412 -0.35355 -0.28336 -0.2939 0.28867 -0.28868
0.16369 0.11363 0.45858 -0.023661 -0.49524 -0.064615 -0.40825 5.9968e-16 0.0015339 -0.063314 -0.49597 0.28868
-0.45418 -0.20003 -0.060894 0.47816 -0.0037812 -0.14611 -0.20412 0.35355 -0.28183 -0.35721 -0.2073 -0.28868
-0.29049 -0.086404 0.39768 0.14424 -0.068707 0.47379 0.20412 0.35355 0.40275 -0.12399 0.2691 0.28868
-0.16369 -0.11363 -0.45858 -0.023661 -0.49524 -0.064615 -0.40825 1.5558e-16 -0.0015339 0.063314 0.49597 0.28868####GULP
0.080147 -0.268189 0.188042 0.188701 -0.376743 -0.188042
-0.188701 0.268189 0.456890 -0.080147 0.376743 -0.456890
-0.036137 0.380558 -0.344420 0.463013 -0.118593 0.344420
-0.463013 -0.380558 0.082456 0.036137 0.118593 -0.082456
-0.492210 0.182347 0.309862 -0.003267 -0.306595 -0.309862
0.003267 -0.182347 -0.185615 0.492210 0.306595 0.185615
-0.029092 0.033400 -0.498034 -0.029092 0.033400 -0.498034
-0.029092 0.033400 -0.498034 -0.029092 0.033400 -0.498034
0.447081 -0.220108 -0.040877 0.447081 -0.220108 -0.040877
0.447081 -0.220108 -0.040877 0.447081 -0.220108 -0.040877
0.221973 0.447702 0.017058 0.221973 0.447702 0.017058
0.221973 0.447702 0.017058 0.221973 0.447702 0.017058
-0.353553 0.353553 -0.000000 0.353553 -0.353553 -0.000000
0.353553 0.353553 0.000000 -0.353553 -0.353553 0.000000
-0.204124 -0.204124 0.408248 0.204124 0.204124 0.408248
0.204124 -0.204124 -0.408248 -0.204124 0.204124 -0.408248
0.231515 0.451818 0.402551 -0.049267 0.171036 -0.402551
0.049267 -0.451818 0.220303 -0.231515 -0.171036 -0.220303
0.363723 -0.081714 0.223299 0.305013 -0.140424 -0.223299
-0.305013 0.081714 -0.445437 -0.363723 0.140424 0.445437
-0.253193 -0.197948 0.195168 0.393115 0.448361 -0.195168
-0.393115 0.197948 0.055246 0.253193 -0.448361 -0.055246
-0.288675 -0.288675 -0.288675 0.288675 0.288675 -0.288675
0.288675 -0.288675 0.288675 -0.288675 0.288675 0.288675The below rij{x,y,z} and Sij{x,y,z} may be helpful
#####rijx
0 0.7818 0.7818 0
-0.7818 0 0 -0.7818
-0.7818 0 0 -0.7818
0 0.7818 0.7818 0#####Rijx
0 0 0 0.7818 0.7818 0.7818 0.7818 0.7818 0.7818 0 0 0
0 0 0 0.7818 0.7818 0.7818 0.7818 0.7818 0.7818 0 0 0
0 0 0 0.7818 0.7818 0.7818 0.7818 0.7818 0.7818 0 0 0
-0.7818 -0.7818 -0.7818 0 0 0 0 0 0 -0.7818 -0.7818 -0.7818
-0.7818 -0.7818 -0.7818 0 0 0 0 0 0 -0.7818 -0.7818 -0.7818
-0.7818 -0.7818 -0.7818 0 0 0 0 0 0 -0.7818 -0.7818 -0.7818
-0.7818 -0.7818 -0.7818 0 0 0 0 0 0 -0.7818 -0.7818 -0.7818
-0.7818 -0.7818 -0.7818 0 0 0 0 0 0 -0.7818 -0.7818 -0.7818
-0.7818 -0.7818 -0.7818 0 0 0 0 0 0 -0.7818 -0.7818 -0.7818
0 0 0 0.7818 0.7818 0.7818 0.7818 0.7818 0.7818 0 0 0
0 0 0 0.7818 0.7818 0.7818 0.7818 0.7818 0.7818 0 0 0
0 0 0 0.7818 0.7818 0.7818 0.7818 0.7818 0.7818 0 0 0###LD.Dynam
75.7050 39.3821 39.3821 -38.3624 -39.3821 0 -38.3624 0 -39.3821 1.0198 0 0
39.3821 75.7050 39.3821 -39.3821 -38.3624 0 0 1.0198 0 0 -38.3624 -39.3821
39.3821 39.3821 75.7050 0 0 1.0198 -39.3821 0 -38.3624 0 -39.3821 -38.3624
-38.3624 -39.3821 0 75.7050 39.3821 -39.3821 1.0198 0 0 -38.3624 0 39.3821
-39.3821 -38.3624 0 39.3821 75.7050 -39.3821 0 -38.3624 39.3821 0 1.0198 0
0 0 1.0198 -39.3821 -39.3821 75.7050 0 39.3821 -38.3624 39.3821 0 -38.3624
-38.3624 0 -39.3821 1.0198 0 0 75.7050 -39.3821 39.3821 -38.3624 39.3821 0
0 1.0198 0 0 -38.3624 39.3821 -39.3821 75.7050 -39.3821 39.3821 -38.3624 0
-39.3821 0 -38.3624 0 39.3821 -38.3624 39.3821 -39.3821 75.7050 0 0 1.0198
1.0198 0 0 -38.3624 0 39.3821 -38.3624 39.3821 0 75.7050 -39.3821 -39.3821
0 -38.3624 -39.3821 0 1.0198 0 39.3821 -38.3624 0 -39.3821 75.7050 39.3821
0 -39.3821 -38.3624 39.3821 0 -38.3624 0 0 1.0198 -39.3821 39.3821 75.7050#####LD.freq
0 + 2.0197i
0 + 2.0197i
0 + 2.0197i
0 + 0.0000i
0.0000
0.0000
8.6421
8.6421
12.3875
12.3875
12.3875
17.6346 #####vijx
1.0e+04 *
Columns 1 through 9
-0.0000 -0.0000 -0.0000 -0.0025 0.1094 - 0.1094i 0.0285 - 0.0285i -0.0000 + 0.0000i -0.0001 + 0.0001i -0.0001 + 0.0001i
-0.0000 0.0000 -0.0000 0.0039 -0.0698 + 0.0698i 0.0484 - 0.0484i -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0002 + 0.0002i
0.0000 0.0000 -0.0000 0.0001 0.0159 - 0.0159i 0.0171 - 0.0171i 0.0001 - 0.0001i 0.0002 - 0.0002i -0.0001 + 0.0001i
0.0025 -0.0039 -0.0001 0.0000 -0.0000 + 0.0000i 0.0000 - 0.0000i 0.5886 - 0.5886i -0.3398 + 0.3398i 0.0316 - 0.0316i
-0.1094 - 0.1094i 0.0698 + 0.0698i -0.0159 - 0.0159i 0.0000 + 0.0000i -0.0000 0.0000 0.0724 -0.0418 -2.5404
-0.0285 - 0.0285i -0.0484 - 0.0484i -0.0171 - 0.0171i -0.0000 - 0.0000i -0.0000 0 -0.0156 0.0090 -0.5210
0.0000 + 0.0000i 0.0000 + 0.0000i -0.0001 - 0.0001i -0.5886 - 0.5886i -0.0724 0.0156 -0.0000 0.0000 -0.0000
0.0001 + 0.0001i 0.0000 + 0.0000i -0.0002 - 0.0002i 0.3398 + 0.3398i 0.0418 -0.0090 -0.0000 -0.0000 0.0000
0.0001 + 0.0001i 0.0002 + 0.0002i 0.0001 + 0.0001i -0.0316 - 0.0316i 2.5404 0.5210 0.0000 -0.0000 -0.0000
0.0002 + 0.0002i -0.0001 - 0.0001i 0.0000 + 0.0000i 0.0350 + 0.0350i -1.1532 0.9743 0.0000 -0.0000 -0.0000
0.0001 + 0.0001i -0.0000 - 0.0000i 0.0000 + 0.0000i 0.0155 + 0.0155i -0.3240 0.6074 -0.0001 0.0000 -0.0000
0.0000 + 0.0000i 0.0000 + 0.0000i -0.0000 - 0.0000i 1.4009 + 1.4009i 0.1723 -0.0372 -0.0000 0.0000 0.0000
Columns 10 through 12
-0.0002 + 0.0002i -0.0001 + 0.0001i -0.0000 + 0.0000i
0.0001 - 0.0001i 0.0000 - 0.0000i -0.0000 + 0.0000i
-0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 - 0.0000i
-0.0350 + 0.0350i -0.0155 + 0.0155i -1.4009 + 1.4009i
1.1532 0.3240 -0.1723
-0.9743 -0.6074 0.0372
-0.0000 0.0001 0.0000
0.0000 -0.0000 -0.0000
0 0.0000 -0.0000
-0.0000 0.0000 -0.0001
-0.0000 0.0000 0.0001
0.0001 -0.0001 -0.0000 #####Sijx
1.0e+05 *
Columns 1 through 9
0 0 - 0.0000i 0 - 0.0000i 0 - 0.0003i 0.0110 + 0.0110i 0.0029 + 0.0029i -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i
0 + 0.0000i 0 0 - 0.0000i 0 + 0.0004i -0.0071 - 0.0071i 0.0049 + 0.0049i -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0001 + 0.0001i
0 - 0.0000i 0 - 0.0000i 0 0 + 0.0000i 0.0016 + 0.0016i 0.0017 + 0.0017i 0.0001 - 0.0000i 0.0001 - 0.0001i -0.0000 + 0.0000i
0 - 0.0003i 0 + 0.0004i 0 + 0.0000i 0 -0.0000 - 0.0000i 0.0000 - 0.0000i 0.2543 - 0.2543i -0.1468 + 0.1468i 0.0195 - 0.0195i
-0.0110 + 0.0110i 0.0071 - 0.0071i -0.0016 + 0.0016i 0.0000 - 0.0000i 0 0.0000 0.0313 -0.0181 -1.5734
-0.0029 + 0.0029i -0.0049 + 0.0049i -0.0017 + 0.0017i -0.0000 - 0.0000i -0.0000 0 -0.0067 0.0039 -0.3227
0.0000 + 0.0000i 0.0000 + 0.0000i -0.0001 - 0.0000i -0.2543 - 0.2543i -0.0313 0.0067 0 0.0000 -0.0000
0.0000 + 0.0000i 0.0000 + 0.0000i -0.0001 - 0.0001i 0.1468 + 0.1468i 0.0181 -0.0039 -0.0000 0 0.0000
0.0000 + 0.0000i 0.0001 + 0.0001i 0.0000 + 0.0000i -0.0195 - 0.0195i 1.5734 0.3227 0.0000 -0.0000 0
0.0001 + 0.0001i -0.0001 - 0.0000i 0.0000 + 0.0000i 0.0217 + 0.0217i -0.7143 0.6034 0.0000 -0.0000 -0.0000
0.0001 + 0.0001i -0.0000 - 0.0000i 0.0000 + 0.0000i 0.0096 + 0.0096i -0.2007 0.3762 -0.0001 0.0000 -0.0000
0.0000 + 0.0000i 0.0000 + 0.0000i -0.0000 - 0.0000i 1.2352 + 1.2352i 0.1520 -0.0328 -0.0000 0.0000 0.0000
Columns 10 through 12
-0.0001 + 0.0001i -0.0001 + 0.0001i -0.0000 + 0.0000i
0.0001 - 0.0000i 0.0000 - 0.0000i -0.0000 + 0.0000i
-0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 - 0.0000i
-0.0217 + 0.0217i -0.0096 + 0.0096i -1.2352 + 1.2352i
0.7143 0.2007 -0.1520
-0.6034 -0.3762 0.0328
-0.0000 0.0001 0.0000
0.0000 -0.0000 -0.0000
0 0.0000 -0.0000
0 0.0000 -0.0001
-0.0000 0 0.0002
0.0001 -0.0002 0 #####Sijx*LD.lor
ans =
Columns 1 through 9
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 0.0000 + 0.0000i 0.0000 + 0.0000i 0 0 0
0 0 0 0.0000 - 0.0000i 0 0.0000 0 0 0
0 0 0 0.0000 - 0.0000i 0.0000 0 0 0 0
0 0 0 0 0 0 0 0.0000 0.0222
0 0 0 0 0 0 0.0000 0 0.0074
0 0 0 0 0 0 0.0222 0.0074 0
0 0 0 0 0 0 0.4396 0.1465 0.0000
0 0 0 0 0 0 1.4194 0.4731 0.0000
0 0 0 0 0 0 0 0 0.1305
Columns 10 through 12
0 0 0
0 0 0
0 0 0
0 0 0
0 0 0
0 0 0
0.4396 1.4194 0
0.1465 0.4731 0
0 0.0000 0.1305
0 0.0000 2.5825
0.0000 0 8.3380
2.5825 8.3380 0 #####Di
ans =
0
0
0
-0.0000 - 0.0000i
0.0000 - 0.0000i
0.0000 - 0.0000i
0.0528
0.0528
0.0925
0.0925
0.0925
0.1116