gulp af

###links http://fortran90.org/src/best-practices.html#modules-and-programs

http://www.csee.umbc.edu/~squire/fortranclass/summary.shtml#Stru

http://fortran90.org/

#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 What are the units of "freq" internally? ###question

####answer The units of freq are in wave numbers (cm^-1) on return from pdiag(g). ##Q 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 And also an option to output the D_i:

output Di Di.gout ###question

##Q

###question

####answer

#result

##GULP

##correct