text:
The UC Davis School of Medicine offers a traditional four year http://www.ucdmc.ucdavis.edu/mdprogram/ M.D. program as well as http://www.ucdmc.ucdavis.edu/medschool/education/ a number of other joint degree and both pre and postgraduate programs.  While the school follows the traditional model of the first two years being class based and the second two years being clinically based, UCDSOM strives to supply its first and second year students with numerous clinical opportunities ranging from http://shifaclinic.org/interclinic/index.html student run clinics to a doctoring curriculum starting in the first quarter of first year.

The School of Medicine moved to a http://www.ucdmc.ucdavis.edu/mdprogram/facilities/index.html new building on the UC Davis Medical Center UCDMC campus in Sacramento in January of 2007. Before the entire program moved to to the Sacramento Medical Center campus, the coursework years of the MD program took place in Davis.

The current dean is Dr. Claire Pomeroy.

 Students

There are approximately 93 med students in each incoming class.  When not found in class or at the hospital, students can be found at student and community run clinics, taking a break from classes at the ARC, playing Ultimate Frisbee, participating in a rowdy game of pingpong in the student lounge, studying, eating at Scrubs or La Bou (or more likely at home), catching some shuteye, and procrastinating by making Front Page Davis Wiki pages.

The atmosphere at UCDSOM tends to be challenging but laid back, most notably during first and second year, since classes are graded on a Pass/Fail basis.  This means better friendships and a noncompetitive, cooperative environment when it comes to school.

 Departments

  Anesthesiology and Pain Medicine
  Biochemistry and Molecular Medicine
  Cell Biology and Human Anatomy
  Dermatology
  Emergency Medicine
  Family and Community Medicine
  Internal Medicine
  Medical Microbiology and Immunology
  Medical Pharmacology and Toxicology
  Neurological Surgery
  Neurology
  Obstetrics and Gynecology
  Ophthalmology and Vision Science
  Orthopaedic Surgery
  Otolaryngology
  Pathology and Laboratory Medicine
  Pharmacology and Toxicology
  Physical Medicine and Rehabilitation
  Physiology and Membrane Biology
  Psychiatry and Behavioral Medicine
  Public Health Sciences
  Radiation Oncology
  Radiology
  Surgery
  Urology

 History

The University of California Regents approved the creation of a school of medicine at UC Davis in 1963.  The school opened its doors to students in September of 1968.  There were 48 students admitted in the first class.  Under the supervision of the founding dean, C. John Tupper, M.D., after whom Tupper Hall is named, it took only 34 months to establish the medical school.  This was and remains the shortest startup time of any school of medicine http://www.ucdmc.ucdavis.edu/pmr/pmr/depthome/1965.html (citation).

The relationship between the school of medicine and the UC Davis Medical Center began in the fall of 1966.  At the time, the facility was the Sacramento County Hospital.  Later it was purchased by the university and assumed its new name.

The School of Medicine was part of the controversy surrounding the historic Supreme Court case of Regents of the University of California v. Bakke, that barred quota systems in affirmative action.  Allan Bakke applied to the School of Medicine and was denied admission.  He sued on the basis that the quotas used in admission were discriminatoryhttp://en.wikipedia.org/wiki/Regents_of_University_of_California_vs._Bakke (citation).

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
        subroutine AsIq (y,       x,       shp,
     &                   shgl,    ien,     xmudmi,
     &                   qres,    rmass    )
c
c----------------------------------------------------------------------
c
c This routine computes and assembles the data corresponding to the
c interior elements for the global reconstruction of the diffusive
c flux vector.
c
c input:
c     y     (nshg,ndof)        : Y variables
c     x     (numnp,nsd)         : nodal coordinates
c     shp   (nshape,ngauss)     : element shape-functions
c     shgl  (nsd,nshape,ngauss) : element local shape-function gradients
c     ien   (npro)              : nodal connectivity array
c
c output:
c     qres  (nshg,nflow-1,nsd)  : residual vector for diffusive flux
c     rmass  (nshg)            : lumped mass matrix
c
c----------------------------------------------------------------------
c
        include "common.h"
c
        dimension y(nshg,ndof),               x(numnp,nsd),
     &            shp(nshl,ngauss),
     &            shgl(nsd,nshl,ngauss),
     &            ien(npro,nshl),
     &            qres(nshg,idflx),    rmass(nshg)
c
c.... element level declarations
c
        dimension ycl(npro,nshl,ndof),        xl(npro,nenl,nsd),
     &            ql(npro,nshl,idflx),       rmassl(npro,nshl),
     &            xmudmi(npro,ngauss)
c
        dimension sgn(npro,nshape)
c
c.... create the matrix of mode signs for the hierarchic basis
c     functions.
c
        if (ipord .gt. 1) then
           call getsgn(ien,sgn)
        endif
c
c.... gather the variables
c

        call localy(y,      ycl,     ien,    ndof,   'gather  ')
        call localx(x,      xl,     ien,    nsd,    'gather  ')
c
c.... get the element residuals
c
        ql     = zero
        rmassl = zero

        call e3q  (ycl,      shp,    shgl,
     &             xl,       ql,     rmassl,
     &             xmudmi,   sgn )

c
c.... assemble the diffusive flux residual
c
        call local (qres,   ql,  ien,  idflx,'scatter ')
        call local (rmass,  rmassl, ien,  1,      'scatter ')
c
c.... end
c
        return
        end

ext: f


text:
!  This is bootstrap unrestricted CCSD program written by the Spring
!  2015 class of Advanced Quantum Mechanics (CHEM 850) at the
!  University of Kansas: Tal Aharon, Matthew Barclay, Sijin Ren, and
!  Marco Caricato.
!
!  This program is designed to accept previously calculated
!  Hartree-Fock two-electron integrals and apply the coupled-cluster
!  post-SCF method.
!
!  This program will be carried out in a number of steps:
!    (1) Reading in the previously calculated integrals
!    (2) Calculating an initial guess of the excitation Amplitudes
!          defined solely as the corresponding two-electron integrals
!    (3) Iteratively solving for the excitation Amplitudes and
!    incorporating them into the final coupled-cluster energy
!
!  This first block creates a Module for Global variables, allowing us
!  to separate our code into a number of Subroutines that each
!  calculate various Terms of the coupled-cluster Amplitude equation.
!
! Input file:
! Method = choice of method (CCD or CCSD)
! EpsI = occupied orbital energies
! EpsA = virtual orbital energies
! abcd = <ab||cd> integrals (in that index order)
! iabc = <ia||bc> integrals (in that index order)
! iajb = <ia||jb> integrals (in that index order)
! ijab = <ij||ab> integrals (in that index order)
! ijka = <ij||ka> integrals (in that index order)
! ijkl = <ij||kl> integrals (in that index order)
!
! The name of the files is generic. If no file name is specified, the
! default names are: EpsA, EpsI, Iabcd, Iiabc, Iiajb, Iijab, Iijka, Iijkl.
!
Module Global
   Double Precision, Allocatable, Dimension(:,:,:,:) :: abcd, iabc, iajb, ijab, ijka, ijkl, DubAmpl, Dijab, AmplDiff
   Double Precision, Allocatable, Dimension (:,:) :: SingAmpl, Dia
   Double Precision, Allocatable, Dimension (:) :: EpsI, EpsA
End Module Global
!
Program Overall
   Implicit None
   Character(300) :: Method, Junk
   Integer IError
   Open(18, File='Input', Status='old', IOStat = IError)
   If (IError .ne. 0) Then
     Print*, "No Input file found"
     Call Exit()
   End If
   Read(18,*) Junk, Junk, Method
   Close(18)
   If (Method == "CCSD") Then
     Call CCSD
   else If (Method == "CCD") Then
     Call CCD
   else
      Write(6,*) "Method Not Recognized"
   End If
End Program Overall
!
!!This next block houses the actual program, which Calls the various Subroutines and calculates the Amplitudes
!    and Energy.
Subroutine CCSD
   Use Global
   Implicit Integer (a-n)
   Character(300) :: FileName, Junk
   Double Precision, Allocatable, Dimension(:,:,:,:) :: UpdDubAmpl
   Double Precision, Allocatable, Dimension(:,:) :: UpdSingAmpl
   Double Precision ConvA, ConvE, Energy, PrevEnergy, Term1, Term2, Term3, Term4, Term5, Term6,&
&Term7, Term8, Term9, Term10, Retval, Ret, Temp, Temp2, Term11, Term12, Term13, Term14,&
&TermA, TermB, TermC, TermD, TermE, TermF, Temp3
   Integer NO, NV, aaa, Counter
   Common NO, NV
!
!!First, the program reads the previously calculated values based on occupied (i,j,k,l...m,n) orbitals and virtual (a,b,c,d...e,f) orbitals:
   Open(40, File='Input')
   Read(40,*)
   Read(40,*,IOStat=IError3) Junk, Junk, FileName
!EpsI  !default if it doesn't find the file is EpsI
   Open(18, File=FileName, Status='old', IOStat=IError)
   If (IError .eq. 0) Then
     Read(18,*) NO
     Allocate(EpsI(NO))
     Do i = 1, NO
        Read(18,*) EpsI(i)
     End Do
     IError2 = 0
   End If
   Close(18)
   If (IError .ne. 0) Then
     Open(18, File='EpsI', IOStat=IError2)
     If (IError2 .eq. 0) Then
       Read(18,*) NO
       Allocate(EpsI(NO))
       Do i = 1, NO
          Read(18,*) EpsI(i)
       End Do
     End If
     Close(18)
   End If
   If ((IError .ne. 0) .and. (IError2 .ne. 0)) Then
     Print*, "EpsI not found"
     Return
   End If
!!!!!where we have Read the orbital energies (EpsI) of the Occupied orbitals
   Read(40,*,IOStat=IError3) Junk, Junk, FileName
   IError2 = 400
!EpsA  !default if it doesn't find the file is EpsA
   Open(9, File=FileName, Status='old', IOStat=IError)
   If (IError .eq. 0) Then
     Read(9,*) NV
     Allocate(EpsA(NV))
     Do i = 1, NV
        Read(9,*) EpsA(i)
     End Do
     IError2=0
   End If
   Close(9)
   If (IError .ne. 0) Then
     Open(9, File='EpsA', IOStat=IError2)
     If (IError2 .eq. 0) Then
       Read(9,*) NV
       Allocate(EpsA(NV))
       Do i = 1, NV
          Read(9,*) EpsA(i)
       End Do
     End If
     Close(9)
   End If
   If ((IError .ne. 0) .and. (IError2 .ne. 0)) Then
     Print*, "EpsA not found"
     Return
   End If
!!!!!similarly, we have Read the energies of the Virtual orbitals
!!!!!Next, the two-electron integrals (of form <wx||yz>) are Read:
   Read(40,*, IOStat=IError3) Junk, Junk, FileName
   IError2 = 400
!Iabcd  !default if it doesn't find the file is Iabcd
   Open(11, File=FileName, Status='old', IOStat=IError)
   If (IError .eq. 0) Then
     Read(11,*) NO, NV
     Allocate(abcd(NV,NV,NV,NV))
     Do d = 1, NV
        Do c = 1, NV
           Do b = 1, NV
              Do a = 1, NV
                 Read(11,*) abcd(a,b,c,d)
              End Do
           End Do
        End Do
     End Do
     IError2 = 0
   End If
   Close(11)
   If (IError .ne. 0) Then
     Open(11, File='Iabcd', Status='old', IOStat=IError2)
     If (IError2 .eq. 0) Then
       Read(11,*) NO, NV
       Allocate(abcd(NV,NV,NV,NV))
       Do d = 1, NV
          Do c = 1, NV
             Do b = 1, NV
                Do a = 1, NV
                   Read(11,*) abcd(a,b,c,d)
                End Do
             End Do
          End Do
       End Do
     End If
     Close(11)
   End If
   If ((IError .ne. 0) .and. (IError2 .ne. 0)) Then
     Print*, "<ab||cd> not found"
     Return
   End If
!!!!!here we have constructed the Global 4-Dimensional matrix of integrals whose indices run over
!         only virtual orbitals (e.g. a,b,c,d)
   Read(40,*, IOStat=IError3) Junk, Junk, FileName
   IError2 = 400
!Iiabc !default if it doesn't find the file is Iiabc
   Open(12, File=FileName, Status='old', IOStat=IError)
   If (IError .eq. 0) Then
     Read(12,*) NO, NV
     Allocate(iabc(NO,NV,NV,NV))
     Do c = 1, NV
        Do b = 1, NV
           Do a = 1, NV
              Do i = 1, NO
                 Read(12,*) iabc(i,a,b,c)
              End Do
           End Do
        End Do
     End Do
     IError2 = 0
   End If
   Close(12)
   If (IError .ne. 0) Then
     Open(12, File='Iiabc', Status='old', IOStat=IError2)
     If (IError2 .eq. 0) Then
       Read(12,*) NO, NV
       Allocate(iabc(NO,NV,NV,NV))
       Do c = 1, NV
          Do b = 1, NV
             Do a = 1, NV
                Do i = 1, NO
                   Read(12,*) iabc(i,a,b,c)
                End Do
             End Do
          End Do
       End Do
     End If
     Close(12)
   End If
   If ((IError .ne. 0) .and. (IError2 .ne. 0)) Then
     Print*, "<ia||bc> not found"
     Return
   End If
!!!!!here we have constructed a similar matrix of integrals involving three virtual orbitals and one occupied orbital
!
!!!!!similarly, the Global variable of two occupied, two virtual integrals in the form <ia||jb>:
   Read(40,*, IOStat=IError3) Junk, Junk, FileName
   IError2 = 400
!Iiajb !default if it doesn't find the file is Iiajb
   Open(13, File=FileName, Status='old', IOStat=IError)
   If (IError .eq. 0) Then
     Read(13,*) NO, NV
     Allocate(iajb(NO,NV,NO,NV))
     Do b = 1, NV
        Do j = 1, NO
           Do a = 1, NV
              Do i = 1, NO
                 Read(13,*) iajb(i,a,j,b)
              End Do
           End Do
        End Do
     End Do
     IError2 = 0
   End If
   Close(13)
   If (IError .ne. 0) Then
     Open(13, File='Iiajb', Status='old', IOStat=IError2)
     If (IError2 .eq. 0) Then
       Read(13,*) NO, NV
       Allocate(iajb(NO,NV,NO,NV))
       Do b = 1, NV
          Do j = 1, NO
             Do a = 1, NV
                Do i = 1, NO
                   Read(13,*) iajb(i,a,j,b)
                End Do
             End Do
          End Do
       End Do
     End If
     Close(13)
   End If
   If ((IError .ne. 0) .and. (IError2 .ne. 0)) Then
     Print*, "<ia||jb> not found"
     Return
   End If
!!!!!and the same for the form <ij||ab>:
   Read(40,*, IOStat=IError3) Junk, Junk, FileName
   IError2 = 400
!Iijab !default if it doesn't find the file is Iijab
   Open(14, File=FileName, Status='old', IOStat=IError)
   If (IError .eq. 0) Then
     Read(14,*) NO, NV
     Allocate(ijab(NO,NO,NV,NV))
     Do b = 1, NV
        Do a = 1, NV
           Do j = 1, NO
              Do i = 1, NO
                 Read(14,*) ijab(i,j,a,b)
              End Do
           End Do
        End Do
     End Do
     IError2 = 0
   End If
   Close(14)
   If (IError .ne. 0) Then
     Open(14, File='Iijab', Status='old', IOStat=IError2)
     If (IError2 .eq. 0) Then
       Read(14,*) NO, NV
       Allocate(ijab(NO,NO,NV,NV))
       Do b = 1, NV
          Do a = 1, NV
             Do j = 1, NO
                Do i = 1, NO
                   Read(14,*) ijab(i,j,a,b)
                End Do
             End Do
          End Do
       End Do
     End If
     Close(14)
   End If
   If ((IError .ne. 0) .and. (IError2 .ne. 0)) Then
     Print*, "<ij||ab> not found"
     Return
   End If
!!!!!and for three occupied, one virtual (<ij||ka>):
   Read(40,*, IOStat=IError3) Junk, Junk, FileName
   IError2 = 400
!Iijka !default if it doesn't find the file is Iijka
   Open(15, File=FileName, Status='old', IOStat=IError)
   If (IError .eq. 0) Then
     Read(15,*) NO, NV
     Allocate(ijka(NO,NO,NO,NV))
     Do a = 1, NV
        Do k = 1, NO
           Do j = 1, NO
              Do i = 1, NO
                 Read(15,*) ijka(i,j,k,a)
              End Do
           End Do
        End Do
     End Do
     IError2 = 0
   End If
   Close(15)
   If (IError .ne. 0) Then
     Open(15, File='Iijka', Status='old', IOStat=IError2)
     If (IError2 .eq. 0) Then
       Read(15,*) NO, NV
       Allocate(ijka(NO,NO,NO,NV))
       Do a = 1, NV
          Do k = 1, NO
             Do j = 1, NO
                Do i = 1, NO
                   Read(15,*) ijka(i,j,k,a)
                End Do
             End Do
          End Do
       End Do
     End If
     Close(15)
   End If
   If ((IError .ne. 0) .and. (IError2 .ne. 0)) Then
     Print*, "<ij||ka> not found"
     Return
   End If
!!!!!and for integrals of only occupied orbitals:
   Read(40,*, IOStat=IError3) Junk, Junk, FileName
   IError2 = 400
!Iijkl !default if it doesn't find the file is Iijkl
   Open(16, File=FileName, Status='old', IOStat=IError)
   If (IError .eq. 0) Then
     Read(16,*) NO, NV
     Allocate(ijkl(NO,NO,NO,NO))
     Do l = 1, NO
        Do k = 1, NO
           Do j = 1, NO
              Do i = 1, NO
                 Read(16,*) ijkl(i,j,k,l)
              End Do
           End Do
        End Do
     End Do
     IError2 = 0
   End If
   Close(16)
   If (IError .ne. 0) Then
     Open(16, File='Iijkl', Status='old', IOStat=IError2)
     If (IError2 .eq. 0) Then
       Read(16,*) NO, NV
       Allocate(ijkl(NO,NO,NO,NO))
       Do l = 1, NO
          Do k = 1, NO
             Do j = 1, NO
                Do i = 1, NO
                   Read(16,*) ijkl(i,j,k,l)
                End Do
             End Do
          End Do
       End Do
     End If
     Close(16)
   End If
   If ((IError .ne. 0) .and. (IError2 .ne. 0)) Then
     Print*, "<ij||kl> not found"
     Return
   End If
   Close(40)
!
!!The next block of the main program constructs the Amplitudes and
!!ultimately deTermines the Energy.
!!!!!First we Allocate memory to the 4-Dimensional matrices for
!!!!!Amplitudes and the denominator Terms 'Dijab' and 'Dia'
!         which (we see below) is the sum of occupied orbital energies
!         minus the virtual orbitals they will be excited to in the
!         corresponding cluster excitation
   Allocate(UpdDubAmpl(NO,NO,NV,NV))
   Allocate(DubAmpl(NO,NO,NV,NV))
   Allocate(Dijab(NO,NO,NV,NV))
   Allocate(AmplDiff(NO,NO,NV,NV))
   Allocate(Dia(NO,NV))
   Allocate(SingAmpl(NO,NV))
   Allocate(UpdSingAmpl(NO,NV))
!
!!!!!Thus we calculate the Dijab and Dia matrices:
   Do b = 1, NV
      Do a = 1, NV
         Do j = 1, NO
            Do i = 1, NO
               Dijab(i,j,a,b) = EpsI(i) + EpsI(j) - EpsA(a) - EpsA(b)
               Dia(i,a) = EpsI(i) - EpsA(a)
            End Do
         End Do
      End Do
   End Do
!
!!Here we perform our initial guess of the excitation Amplitudes based
!!solely on the two-electron integrals from
!    the SCF calculation:
   Do b = 1, NV
      Do a = 1, NV
         Do j = 1, NO
            SingAmpl(j,a) = 0
            Do i = 1, NO
               DubAmpl(i,j,a,b) = ijab(i,j,a,b)/Dijab(i,j,a,b)
            End Do
         End Do
      End Do
   End Do
!!!!!We Then perform an initial Energy calculation based on these
!!!!!Amplitude guesses, giving the MP2 energy:
   Do b = 1, NV
      Do a = 1, NV
         Do j = 1, NO
            Do i = 1, NO
               Call Tau(DubAmpl(i,j,a,b),i,j,a,b, Retval)
               Energy = Energy + .25*Retval*ijab(i,j,a,b)
!!!!!where 'Tau' is a Subroutine function of the excitation Amplitudes
            End Do
         End Do
      End Do
   End Do
   print*, "MP2 Energy", Energy
!
!!The program now proceeds to perform the actual iterations.
!!!!!Each of the numbered "Terms" (Term1, Term2, Term3...) is one of
!!!!!the parts of the Double-excitation cluster operator (T2)
!!!!!equation, as it is broken Down and explicitly defined before
!!!!!being recombined at the End.
   Energy = 0
   PrevEnergy = 0
   aaa = 0
   Counter = 0
!!!!!thus we have provided an initial guess of zero for each of our
!!!!!following variables, allowing for ease of summation
!
!
!!!!!The following block of code is the main structure of the
!!!!!iterative loop,including the test of
!         and criteria for Convergence.
   Do
      If (Counter == 100) Then
         print*, "Convergence not achieved. Current energy after 100 cycles is ", Energy
         Exit
      End If
      If (aaa == 1) Then !Convergence
         Exit
      End If
!!!!!where we have initially determined that If the program Does not Converge within 100 iterations, it will cease.
      Do b = 1, NV
         Do a = 1, NV
            Do j = 1, NO
               Do i = 1, NO
                  Term1 = 0
                  Term2 = 0
                  Term3 = 0
                  Term4 = 0
                  Term5 = 0
                  Term6 = 0
                  Term7 = 0
                  Term8 = 0
                  Term9 = 0
                  Term10 = 0
                  Term11 = 0
                  Term12 = 0
                  Term13 = 0
                  Term14 = 0
!!!!!here we have constructed the iterative loops over all indices of occupied (i,j) and virtual (a,b) orbitals
!         prior to actually defining the individual Terms
!!Terms 1 and 2 consist of two permutations of the virtual orbitals (ab) within a sum over a single virtual orbital
!    (e) of Terms of the "F" matrix [defined as a Subroutine] whose corresponding indices are virtual orbitals:
                  Do e = 1, NV
                     Temp = 0
                     Temp2 = 0
                     Do m = 1, NO
                        Call FBoth(m, e, Retval)
                        Temp = Temp - ((.5)*SingAmpl(m,b)*Retval)
                        Temp2 = Temp2 - ((.5) * SingAmpl(m,a)*Retval)
!mc                        Temp2 = Temp2 + ((.5) * SingAmpl(m,a)*Retval)
                     End Do
                     Call FVirt(b,e,Retval)
                     Term1 = Term1 + DubAmpl(i,j,a,e) * Retval + DubAmpl(i,j,a,e)*Temp
                     Call FVirt(a,e,Retval)
                     Term2 = Term2 + DubAmpl(i,j,b,e) * Retval + DubAmpl(i,j,b,e)*Temp2
                  End Do
!!!!!thus we have summed every (*,e) element of the virtual-orbital F matrix ['FVirt'], constructing the Terms 1 and 2
!         [F(b,e) and F(a,e) respectively] alongside the corresponding excitation Amplitudes
!!Similarly, Terms 3 and 4 consist of the contributions where a single occupied orbital (m) dIffers in Terms containing
!    the occupied-orbital F matrix ['FOcc']; each Term is a permutation of the occupied orbitals (i,j):
                  Do m = 1, NO
                     Temp = 0
                     Temp2 = 0
                     Do e = 1, NV
                        Call FBoth(m, e, Retval)
                        Temp = Temp + ((.5)*SingAmpl(j,e)*Retval)
                        Temp2 = Temp2 + ((.5) * SingAmpl(i,e)*Retval)
!mc                        Temp2 = Temp2 - ((.5) * SingAmpl(i,e)*Retval)
                     End Do
                     Call FOcc(m, j, Retval)
                     Term3 = Term3 + DubAmpl(i,m,a,b) * Retval + DubAmpl(i,m,a,b) * Temp
                     Call FOcc(m, i, Retval)
                     Term4 = Term4 + DubAmpl(j,m,a,b) * Retval + DubAmpl(j,m,a,b) * Temp2
                  End Do
!!Terms 5 and 6 are simple sums over two occupied (m,n) and two virtual (e,f) orbitals respectively, and Call upon the
!    corresponding W matrices and excitation Amplitudes:
                  Do n = 1, NO
                     Do m = 1, NO
                        Call Tau(DubAmpl(m,n,a,b),m,n,a,b, Ret)
                        Call WOcc(m,n,i,j, Retval)
                        Term5 = Term5 + .5*Ret*Retval
                     End Do
                  End Do
                  Do f = 1, NV
                     Do e = 1, NV
                        Call Tau(DubAmpl(i,j,e,f),i,j,e,f, Ret)
                        Call WVirt(a,b,e,f, Retval)
                        Term6 = Term6 + .5*Ret*Retval
                     End Do
                  End Do
!!Terms 6-10 are the contributions of the mixed (two occupied, two virtual) W matrix ['WBoth'] summed over one dIffering
!    occupied and one dIffering virtual orbital; each Term is part of simultaneous permutations of (a,b) and (i,j):
                  Do e = 1, NV
                     Do m = 1, NO
                        Call WBoth(m,b,e,j, Retval)
                        Term7 = Term7 + DubAmpl(i,m,a,e)*Retval - SingAmpl(i,e)*SingAmpl(m,a)&
                        &*(-iajb(m,b,j,e))
                        Call WBoth(m,a,e,j, Retval)
                        Term8 = Term8 + DubAmpl(i,m,b,e)*Retval  - SingAmpl(i,e)*SingAmpl(m,b)&
                        &*(-iajb(m,a,j,e))
                        Call WBoth(m,b,e,i, Retval)
                        Term9 = Term9 + DubAmpl(j,m,a,e)*Retval - SingAmpl(j,e)*SingAmpl(m,a)&
                        &*(-iajb(m,b,i,e))
                        Call WBoth(m,a,e,i, Retval)
                        Term10 = Term10 + DubAmpl(j,m,b,e)*Retval - SingAmpl(j,e)*SingAmpl(m,b)&
                        &*(-iajb(m,a,i,e))
                     End Do
                  End Do
!
                  Do e = 1, NV
                     Term11 = Term11 + SingAmpl(i,e) * iabc(j,e,b,a)
                     Term12 = Term12 + SingAmpl(j,e) * iabc(i,e,b,a)
                  End Do
!
                  Do m = 1, NO
                     Term13 = Term13 + SingAmpl(m,a) * ijka(i,j,m,b)
                     Term14 = Term14 + SingAmpl(m,b) * ijka(i,j,m,a)
                  End Do
!
!!Having defined all the Terms in the Amplitude equation, we combine them all together in the 'updated' excitation
!    Amplitude:
                  UpdDubAmpl(i,j,a,b) = (ijab(i,j,a,b)+Term1-Term2-Term3&
&+Term4+Term5+Term6+Term7-Term8-Term9+Term10+Term11-Term12-Term13+Term14)/Dijab(i,j,a,b)
              End Do
            End Do
         End Do
      End Do
!Here the Terms for the T1 Amplitude are calculated, similar to those of the T2 Amplitudes
      Do a = 1, NV
         Do i = 1, NO
            TermA = 0
            TermB = 0
            TermC = 0
            TermD = 0
            TermE = 0
            TermF = 0
            Do e = 1, NV
               Call FVirt(a,e,Retval)
               TermA = TermA + SingAmpl(i,e)*Retval
            End Do
            Do e = 1, NV
               Do m = 1, NO
                  Call FBoth(m,e,Retval)
                  TermC = TermC + DubAmpl(i,m,a,e)*Retval
               End Do
            End Do
            Do m = 1, NO
               Call FOcc(m,i,Retval)
               TermB = TermB + SingAmpl(m,a) * Retval
            End Do
            Do f = 1, NV
               Do n = 1, NO
                  TermD = TermD + SingAmpl(n,f) * iajb(n,a,i,f)
               End Do
            End Do
            Do f = 1, NV
               Do e = 1, NV
                  Do m = 1, NO
                     TermE = TermE + .5*DubAmpl(i,m,e,f)*iabc(m,a,e,f)
                  End Do
               End Do
            End Do
            Do n = 1, NO
               Do e = 1, NV
                  Do m = 1, NO
                     TermF = TermF + .5*DubAmpl(m,n,a,e)*ijka(m,n,i,e)
                  End Do
               End Do
            End Do
            UpdSingAmpl(i,a) = (TermA - TermB + TermC - TermD - TermE - TermF)/Dia(i,a)
        End Do
      End Do
      AmplDiff = UpdDubAmpl-DubAmpl
      !print*, UpdDubAmpl(5,6,2,3), " ", DubAmpl(5,6,2,3), " ", AmplDiff(5,6,2,3)
      DubAmpl = UpdDubAmpl
      SingAmpl = UpdSingAmpl
!!!!!and thus the Previous excitation Amplitude is redefined as the updated value
!
!!Finally, the program recalculates the coupled-cluster Energy from the new Amplitudes:
      Energy = 0
      Do b = 1, NV
         Do a = 1, NV
            Do j = 1, NO
               Do i = 1, NO
                  Call Tau(DubAmpl(i,j,a,b),i,j,a,b, Retval)
                  Energy = Energy + .25*Retval*ijab(i,j,a,b)
               End Do
            End Do
         End Do
      End Do
!
!!and checks for Convergence based on both the Energy dIfference (between current Energy and Previous Energy), as well
!    as the dIfference in the root-mean squared (RMS) value of the excitation Amplitudes:
      Call AmplConv(ConvA)
      ConvE = Energy - PrevEnergy
      ConvE = sqrt(ConvE**2)
      If ((ConvA <= 1.0D-5) .and. (ConvE <= 1.0D-7)) Then
         aaa = 1
      End If
!!!!!If Convergence has been achieved, the program signals the End of the iterative loop and prints the final Energy;
!         otherwise, it cycles back through the next iteration
      PrevEnergy = Energy
      Counter = Counter + 1
      Print*, "The energy of iteration ", Counter, "is", Energy
   End Do
   print*, "The final Energy is: ", Energy
End Subroutine CCSD
!
!
!
!!The following blocks are each of the Subroutines/Terms incorporated into the overarching Amplitude equation.
!
!
!!The Tau and "TauTilde" Subroutines calculate the values of the variables of their namesake, which are basiCally sums
!    of the operators leading to Double-excitation.
!!!!!As this program's current function is a CCD calculation, the current values of Tau/Tautilde are solely that of the
!         corresponding Double-excitation Amplitude.
Subroutine Tau(Tijab, i, j, a, b, Retval)
   Use Global
   Implicit None
   Double Precision, Intent(in) :: Tijab
   Integer, Intent(in) :: i, a, j, b
   Double Precision, Intent(out) :: Retval
   Double Precision SingAmpDIff
   SingAmpDIff = SingAmpl(i,a)*SingAmpl(j,b) - SingAmpl(i,b)*SingAmpl(j,a)
   Retval = Tijab + SingAmpDIff
End Subroutine Tau
!
Subroutine TauTilde(Tijab, i, j, a, b, Retval)
   Use Global
   Implicit None
   Double Precision, Intent(in) :: Tijab
   Integer, Intent(in) :: i, a, j, b
   Double Precision, Intent(out) :: Retval
   Double Precision SingAmpDIff
   SingAmpDIff = (SingAmpl(i,a)*SingAmpl(j,b) - SingAmpl(i,b)*SingAmpl(j,a))/2.0
   Retval = Tijab + SingAmpDIff
End Subroutine TauTilde
!
!!The F matrices similarly correspond to their namesake and are defined by the nature of their indices - e.g. 'FVirt'
!         is the F matrix defined only through the virtual orbitals of the excitation.
Subroutine FVirt(b, e, Retval)
   Use Global
   Implicit Integer (f-o)
   Integer, Intent(in) :: b, e
   Double Precision, Intent(out) :: Retval
   Double Precision SumMNF, TaoT, SumMF
   Integer NO, NV
   Common NO, NV
   SumMNF = 0
   SumMF = 0
   Do f = 1, NV
      Do m = 1, NO
         SumMF = SumMF + SingAmpl(m,f) * iabc(m,b,f,e)
      End Do
   End Do
   Do f = 1, NV
      Do n = 1, NO
         Do m = 1, NO
            Call TauTilde(DubAmpl(m,n,b,f),m,n,b,f, TaoT)
            SumMNF = SumMNF + TaoT*ijab(m,n,e,f)
         End Do
      End Do
   End Do
   Retval = SumMNF * (-.5) + SumMF
End Subroutine FVirt
!
Subroutine FOcc(m, j, Retval)
   Use Global
   Implicit None
   Integer, Intent(in) :: j, m
   Double Precision, Intent(out) :: Retval
   Double Precision SumNEF, TaoT, SumNE
   Integer n, e, f, NO, NV
   Common NO, NV
   SumNE = 0
   SumNEF = 0
   Do e = 1, NV
      Do n = 1, NO
         SumNE = SumNE + SingAmpl(n,e) * ijka(m,n,j,e)
      End Do
   End Do
   Do f = 1, NV
      Do e = 1, NV
         Do n = 1, NO
            Call TauTilde(DubAmpl(j,n,e,f),j,n,e,f, TaoT)
            SumNEF = SumNEF + TaoT*ijab(m,n,e,f)
         End Do
      End Do
   End Do
   Retval = SumNEF * .5 + SumNE
End Subroutine FOcc
!
Subroutine FBoth(m, e, Retval)
   Use Global
   Implicit None
   Integer, Intent(in) :: e, m
   Double Precision, Intent(out) :: Retval
   Double Precision SumNF
   Integer n, f, NO, NV
   Common NO, NV
   SumNF = 0
   Do f = 1, NV
      Do n = 1, NO
         SumNF = SumNF + SingAmpl(n,f) * ijab(m,n,e,f)
      End Do
   End Do
   Retval = SumNF
End Subroutine FBoth
!
!
!!The W matrices similarly correspond to their namesake and are defined by the nature of their indices - e.g. 'WVirt'
!         is the W matrix defined only through the virtual orbitals of the excitation.
Subroutine WOcc(m,n,i,j,Retval)
   Use Global
   Implicit None
   Integer, Intent(in) :: m, n, i, j
   Double Precision, Intent(out) :: Retval
   Double Precision SumEF, Tao, SumE, SumEPerm
   Integer e, f, NO, NV
   Common NO, NV
   SumEF = 0
   SumE = 0
   SumEPerm = 0
   Do e = 1, NV
      SumE = SumE + SingAmpl(j,e) * (ijka(m,n,i,e))
      SumEPerm = SumEPerm - (SingAmpl(i,e) * (ijka(m,n,j,e)))
   End Do
   Do f = 1, NV
      Do e = 1, NV
         Call Tau(DubAmpl(i,j,e,f),i,j,e,f, Tao)
         SumEF = SumEF +.25*Tao*ijab(m,n,e,f)
      End Do
   End Do
   Retval = SumEF + ijkl(m,n,i,j) + SumE + SumEPerm
End Subroutine WOcc
!
Subroutine WVirt(a,b,e,f,Retval)
   Use Global
   Implicit None
   Integer, Intent(in) :: a, b, e, f
   Double Precision, Intent(out) :: Retval
   Double Precision SumMN, Tao, SumM, SumMPerm
   Integer m, n, NO, NV
   Common NO, NV
   SumMN = 0
   SumM = 0
   SumMPerm = 0
   Do m = 1, NO
      SumM = SumM - SingAmpl(m,b) * (-iabc(m,a,e,f))
      SumMPerm = SumMPerm + (SingAmpl(m,a) * (-iabc(m,b,e,f)))
   End Do
   Do n = 1, NO
      Do m = 1, NO
         Call Tau(DubAmpl(m,n,a,b),m,n,a,b, Tao)
         SumMN = SumMN +.25*Tao*ijab(m,n,e,f)
      End Do
   End Do
   Retval = SumMN + abcd(a,b,e,f) + SumM + SumMPerm
End Subroutine WVirt
!
Subroutine WBoth(m,b,e,j,Retval)
   Use Global
   Implicit None
   Integer, Intent(in) :: m, b, e, j
   Double Precision, Intent(out) :: Retval
   Double Precision SumNF, AmpTerm, SumN, SumF
   Integer n, f, NO, NV
   Common NO, NV
   SumNF = 0
   SumN = 0
   SumF = 0
   Do f = 1, NV
      SumF = SumF + SingAmpl(j,f) * iabc(m,b,e,f)
   End Do
   Do n = 1, NO
      SumN = SumN - (SingAmpl(n,b) * ijka(n,m,j,e))
   End Do
   Do f = 1, NV
      Do n = 1, NO
         AmpTerm = (.5)*DubAmpl(j,n,f,b) + SingAmpl(j,f)*SingAmpl(n,b)
         SumNF = SumNF - AmpTerm*ijab(m,n,e,f)
!mc         AmpTerm = (-.5)*DubAmpl(j,n,f,b) + SingAmpl(j,f)*SingAmpl(n,b)
!mc         SumNF = SumNF + AmpTerm*ijab(m,n,e,f)
      End Do
   End Do
   Retval = SumNF + (-iajb(m,b,j,e)) + SumF + SumN
End Subroutine WBoth
!
!This program is designed to accept previously calculated Hartree-Fock two-electron integrals and apply the
!    coupled-cluster post-SCF method. Currently, the program only incorporates contributions of Double excitations (CCD)
!This program will be carried out in a number of steps:
!    (1) Reading in the previously calculated integrals
!    (2) calculating an initial guess of the excitation amplitudes defined solely as the corresponding
!          two-electron integrals
!    (3) iteratively solving for the excitation amplitudes and incorporating them into the final coupled-cluster Energy
!!This next block hoUses the actual program, which Calls the various Subroutines and calculates the amplitudes
!    and Energy.
!
Subroutine CCD
   Use Global
   Implicit Integer (a-n)
   Double Precision, Allocatable, Dimension(:,:,:,:) :: UpdDubAmpl
   Character(300) :: FileName, Junk
   Double Precision ConvA, ConvE, Energy, PrevEnergy, Term1, Term2, Term3, Term4, Term5, Term6,&
&Term7, Term8, Term9, Term10, Retval, Ret
   Integer NO, NV, aaa, Counter
   Common NO, NV
!!First, the program Reads the previously calculated values based on occupied (i,j,k,l...m,n) orbitals and
!    virtual (a,b,c,d...e,f) orbitals:
   Open(40, File='Input')
   Read(40,*)
   Read(40,*,IOStat=IError3) Junk, Junk, FileName
!EpsI  !default if it doesn't find the file is EpsI
   Open(18, File=FileName, Status='old', IOStat=IError)
   If (IError .eq. 0) Then
     Read(18,*) NO
     Allocate(EpsI(NO))
     Do i = 1, NO
        Read(18,*) EpsI(i)
     End Do
     IError2 = 0
   End If
   Close(18)
   If (IError .ne. 0) Then
     Open(18, File='EpsI', IOStat=IError2)
     If (IError2 .eq. 0) Then
       Read(18,*) NO
       Allocate(EpsI(NO))
       Do i = 1, NO
          Read(18,*) EpsI(i)
       End Do
     End If
     Close(18)
   End If
   If ((IError .ne. 0) .and. (IError2 .ne. 0)) Then
     Print*, "EpsI not found"
     Return
   End If
!!!!!where we have Read the orbital energies (EpsI) of the Occupied orbitals
   Read(40,*,IOStat=IError3) Junk, Junk, FileName
   IError2 = 400
!EpsA  !default if it doesn't find the file is EpsA
   Open(9, File=FileName, Status='old', IOStat=IError)
   If (IError .eq. 0) Then
     Read(9,*) NV
     Allocate(EpsA(NV))
     Do i = 1, NV
        Read(9,*) EpsA(i)
     End Do
     IError2=0
   End If
   Close(9)
   If (IError .ne. 0) Then
     Open(9, File='EpsA', IOStat=IError2)
     If (IError2 .eq. 0) Then
       Read(9,*) NV
       Allocate(EpsA(NV))
       Do i = 1, NV
          Read(9,*) EpsA(i)
       End Do
     End If
     Close(9)
   End If
   If ((IError .ne. 0) .and. (IError2 .ne. 0)) Then
     Print*, "EpsA not found"
     Return
   End If
!!!!!similarly, we have Read the energies of the Virtual orbitals
!!!!!Next, the two-electron integrals (of form <wx||yz>) are Read:
   Read(40,*, IOStat=IError3) Junk, Junk, FileName
   IError2 = 400
!Iabcd  !default if it doesn't find the file is Iabcd
   Open(11, File=FileName, Status='old', IOStat=IError)
   If (IError .eq. 0) Then
     Read(11,*) NO, NV
     Allocate(abcd(NV,NV,NV,NV))
     Do d = 1, NV
        Do c = 1, NV
           Do b = 1, NV
              Do a = 1, NV
                 Read(11,*) abcd(a,b,c,d)
              End Do
           End Do
        End Do
     End Do
     IError2 = 0
   End If
   Close(11)
   If (IError .ne. 0) Then
     Open(11, File='Iabcd', Status='old', IOStat=IError2)
     If (IError2 .eq. 0) Then
       Read(11,*) NO, NV
       Allocate(abcd(NV,NV,NV,NV))
       Do d = 1, NV
          Do c = 1, NV
             Do b = 1, NV
                Do a = 1, NV
                   Read(11,*) abcd(a,b,c,d)
                End Do
             End Do
          End Do
       End Do
     End If
     Close(11)
   End If
   If ((IError .ne. 0) .and. (IError2 .ne. 0)) Then
     Print*, "<ab||cd> not found"
     Return
   End If
!!!!!here we have constructed the Global 4-Dimensional matrix of integrals whose indices run over
!         only virtual orbitals (e.g. a,b,c,d)
   Read(40,*, IOStat=IError3) Junk, Junk, FileName
   IError2 = 400
!Iiabc  !default if it doesn't find the file is Iiabc
   Open(12, File=FileName, Status='old', IOStat=IError)
   If (IError .eq. 0) Then
     Read(12,*) NO, NV
     Allocate(iabc(NO,NV,NV,NV))
     Do c = 1, NV
        Do b = 1, NV
           Do a = 1, NV
              Do i = 1, NO
                 Read(12,*) iabc(i,a,b,c)
              End Do
           End Do
        End Do
     End Do
     IError2 = 0
   End If
   Close(12)
   If (IError .ne. 0) Then
     Open(12, File='Iiabc', Status='old', IOStat=IError2)
     If (IError2 .eq. 0) Then
       Read(12,*) NO, NV
       Allocate(iabc(NO,NV,NV,NV))
       Do c = 1, NV
          Do b = 1, NV
             Do a = 1, NV
                Do i = 1, NO
                   Read(12,*) iabc(i,a,b,c)
                End Do
             End Do
          End Do
       End Do
     End If
     Close(12)
   End If
   If ((IError .ne. 0) .and. (IError2 .ne. 0)) Then
     Print*, "<ia||bc> not found"
     Return
   End If
!!!!!here we have constructed a similar matrix of integrals involving three virtual orbitals and one occupied orbital
!
!!!!!similarly, the Global variable of two occupied, two virtual integrals in the form <ia||jb>:
   Read(40,*, IOStat=IError3) Junk, Junk, FileName
   IError2 = 400
!Iiajb  !default if it doesn't find the file is Iiajb
   Open(13, File=FileName, Status='old', IOStat=IError)
   If (IError .eq. 0) Then
     Read(13,*) NO, NV
     Allocate(iajb(NO,NV,NO,NV))
     Do b = 1, NV
        Do j = 1, NO
           Do a = 1, NV
              Do i = 1, NO
                 Read(13,*) iajb(i,a,j,b)
              End Do
           End Do
        End Do
     End Do
     IError2 = 0
   End If
   Close(13)
   If (IError .ne. 0) Then
     Open(13, File='Iiajb', Status='old', IOStat=IError2)
     If (IError2 .eq. 0) Then
       Read(13,*) NO, NV
       Allocate(iajb(NO,NV,NO,NV))
       Do b = 1, NV
          Do j = 1, NO
             Do a = 1, NV
                Do i = 1, NO
                   Read(13,*) iajb(i,a,j,b)
                End Do
             End Do
          End Do
       End Do
     End If
     Close(13)
   End If
   If ((IError .ne. 0) .and. (IError2 .ne. 0)) Then
     Print*, "<ia||jb> not found"
     Return
   End If
!!!!!and the same for the form <ij||ab>:
   Read(40,*, IOStat=IError3) Junk, Junk, FileName
   IError2 = 400
!Iijab
   Open(14, File=FileName, Status='old', IOStat=IError)
   If (IError .eq. 0) Then
     Read(14,*) NO, NV
     Allocate(ijab(NO,NO,NV,NV))
     Do b = 1, NV
        Do a = 1, NV
           Do j = 1, NO
              Do i = 1, NO
                 Read(14,*) ijab(i,j,a,b)
              End Do
           End Do
        End Do
     End Do
     IError2 = 0
   End If
   Close(14)
   If (IError .ne. 0) Then
     Open(14, File='Iijab', Status='old', IOStat=IError2)
     If (IError2 .eq. 0) Then
       Read(14,*) NO, NV
       Allocate(ijab(NO,NO,NV,NV))
       Do b = 1, NV
          Do a = 1, NV
             Do j = 1, NO
                Do i = 1, NO
                   Read(14,*) ijab(i,j,a,b)
                End Do
             End Do
          End Do
       End Do
     End If
     Close(14)
   End If
   If ((IError .ne. 0) .and. (IError2 .ne. 0)) Then
     Print*, "<ij||ab> not found"
     Return
   End If
!!!!!and for three occupied, one virtual (<ij||ka>):
   Read(40,*, IOStat=IError3) Junk, Junk, FileName
   IError2 = 400
!Iijka  !default if it doesn't find the file is Iijka
   Open(15, File=FileName, Status='old', IOStat=IError)
   If (IError .eq. 0) Then
     Read(15,*) NO, NV
     Allocate(ijka(NO,NO,NO,NV))
     Do a = 1, NV
        Do k = 1, NO
           Do j = 1, NO
              Do i = 1, NO
                 Read(15,*) ijka(i,j,k,a)
              End Do
           End Do
        End Do
     End Do
     IError2 = 0
   End If
   Close(15)
   If (IError .ne. 0) Then
     Open(15, File='Iijka', Status='old', IOStat=IError2)
     If (IError2 .eq. 0) Then
       Read(15,*) NO, NV
       Allocate(ijka(NO,NO,NO,NV))
       Do a = 1, NV
          Do k = 1, NO
             Do j = 1, NO
                Do i = 1, NO
                   Read(15,*) ijka(i,j,k,a)
                End Do
             End Do
          End Do
       End Do
     End If
     Close(15)
   End If
   If ((IError .ne. 0) .and. (IError2 .ne. 0)) Then
     Print*, "<ij||ka> not found"
     Return
   End If
!!!!!and for integrals of only occupied orbitals:
   Read(40,*, IOStat=IError3) Junk, Junk, FileName
   IError2 = 400
!Iijkl  !default if it doesn't find the file is Iijkl
   Open(16, File=FileName, Status='old', IOStat=IError)
   If (IError .eq. 0) Then
     Read(16,*) NO, NV
     Allocate(ijkl(NO,NO,NO,NO))
     Do l = 1, NO
        Do k = 1, NO
           Do j = 1, NO
              Do i = 1, NO
                 Read(16,*) ijkl(i,j,k,l)
              End Do
           End Do
        End Do
     End Do
     IError2 = 0
   End If
   Close(16)
   If (IError .ne. 0) Then
     Open(16, File='Iijkl', Status='old', IOStat=IError2)
     If (IError2 .eq. 0) Then
       Read(16,*) NO, NV
       Allocate(ijkl(NO,NO,NO,NO))
       Do l = 1, NO
          Do k = 1, NO
             Do j = 1, NO
                Do i = 1, NO
                   Read(16,*) ijkl(i,j,k,l)
                End Do
             End Do
          End Do
       End Do
     End If
     Close(16)
   End If
   If ((IError .ne. 0) .and. (IError2 .ne. 0)) Then
     Print*, "<ij||kl> not found"
     Return
   End If
   Close(40)
!!The next block of the main program constructs the amplitudes and ultimately deTermines the Energy.
!!!!!First we Allocate memory to the 4-Dimensional matrices for amplitudes and the denominator Term 'Dijab'
!         which (we see below) is the sum of occupied orbital energies minus the virtual orbitals they will be
!         excited to in the corresponding cluster excitation
   Allocate(UpdDubAmpl(NO,NO,NV,NV))
   Allocate(DubAmpl(NO,NO,NV,NV))
   Allocate(Dijab(NO,NO,NV,NV))
   Allocate(AmplDiff(NO,NO,NV,NV))
!!!!!Thus we calculate the Dijab matrix:
   Do b = 1, NV
      Do a = 1, NV
         Do j = 1, NO
            Do i = 1, NO
               Dijab(i,j,a,b) = EpsI(i) + EpsI(j) - EpsA(a) - EpsA(b)
            End Do
         End Do
      End Do
   End Do
!!Here we perform our initial guess of the excitation amplitudes based solely on the two-electron integrals from
!    the SCF calculation:
   Do b = 1, NV
      Do a = 1, NV
         Do j = 1, NO
            Do i = 1, NO
               DubAmpl(i,j,a,b) = ijab(i,j,a,b)/Dijab(i,j,a,b)
            End Do
         End Do
      End Do
   End Do
!!!!!We then perform an initial Energy calculation based on these amplitude guesses:
   Do b = 1, NV
      Do a = 1, NV
         Do j = 1, NO
            Do i = 1, NO
               Energy = Energy + .25*DubAmpl(i,j,a,b)*ijab(i,j,a,b)
            End Do
         End Do
      End Do
   End Do
   print*, "MP2 Energy", Energy
!!The program now proceeds to perform the actual iterations.
!!!!!Each of the "Terms" is one of the parts of the Double-excitation cluster operator (T2) equation, as it is
!         broken Down and explicitly defined before being recombined at the End.
   Energy = 0
   PrevEnergy = 0
   aaa = 0
   Counter = 0
!!!!!thus we have provided an initial guess of zero for each of our following variables, allowing for ease of summation
!!!!!The following block of code is the main structure of the iterative loop,including the test of
!         and criteria for Convergence.
   Do
      If (Counter == 100) Then
         print*, "Convergence not achieved. Current energy after 100 cycles is ", Energy
         Exit
      End If
      If (aaa == 1) Then !Convergence
         Exit
      End If
!!!!!where we have initially deTermined that If the program Does not Converge within 100 iterations, it will cease.
      Do b = 1, NV
         Do a = 1, NV
            Do j = 1, NO
               Do i = 1, NO
                  Term1 = 0
                  Term2 = 0
                  Term3 = 0
                  Term4 = 0
                  Term5 = 0
                  Term6 = 0
                  Term7 = 0
                  Term8 = 0
                  Term9 = 0
                  Term10 = 0
!!!!!here we have constructed the iterative loops over all indices of occupied (i,j) and virtual (a,b) orbitals
!         prior to actually defining the individual Terms
!!Terms 1 and 2 consist of two permutations of the virtual orbitals (ab) within a sum over a single virtual orbital
!    (e) of Terms of the "F" matrix [defined as a Subroutine] whose corresponding indices are virtual orbitals:
                  Do e = 1, NV
                     Call FV(b,e,Retval)
                     Term1 = Term1 + DubAmpl(i,j,a,e) * Retval
                     Call FV(a,e,Retval)
                     Term2 = Term2 + DubAmpl(i,j,b,e) * Retval
                  End Do
!!!!!thus we have summed every (*,e) element of the virtual-orbital F matrix ['FVirt'], constructing the Terms 1 and 2
!         [F(b,e) and F(a,e) respectively] alongside the corresponding excitation amplitudes
!!Similarly, Terms 3 and 4 consist of the contributions where a single occupied orbital (m) dIffers in Terms containing
!    the occupied-orbital F matrix ['FOcc']; each Term is a permutation of the occupied orbitals (i,j):
                  Do m = 1, NO
                     Call FO(m, j, Retval)
                     Term3 = Term3 + DubAmpl(i,m,a,b) * Retval
                     Call FO(m, i, Retval)
                     Term4 = Term4 + DubAmpl(j,m,a,b) * Retval
                  End Do
!!Terms 5 and 6 are simple sums over two occupied (m,n) and two virtual (e,f) orbitals respectively, and Call upon the
!    corresponding W matrices and excitation amplitudes:
                  Do n = 1, NO
                     Do m = 1, NO
                        Call WO(m,n,i,j, Retval)
                        Term5 = Term5 + .5*DubAmpl(m,n,a,b)*Retval
                     End Do
                  End Do
                  Do f = 1, NV
                     Do e = 1, NV
                        Call WV(a,b,e,f, Retval)
                        Term6 = Term6 + .5*DubAmpl(i,j,e,f)*Retval
                     End Do
                  End Do
!!Terms 6-10 are the contributions of the mixed (two occupied, two virtual) W matrix ['WBoth'] summed over one dIffering
!    occupied and one dIffering virtual orbital; each Term is part of simultaneous permutations of (a,b) and (i,j):
                  Do e = 1, NV
                     Do m = 1, NO
                        Call WB(m,b,e,j, Retval)
                        Term7 = Term7 + DubAmpl(i,m,a,e)*Retval
                        Call WB(m,a,e,j, Retval)
                        Term8 = Term8 + DubAmpl(i,m,b,e)*Retval
                        Call WB(m,b,e,i, Retval)
                        Term9 = Term9 + DubAmpl(j,m,a,e)*Retval
                        Call WB(m,a,e,i, Retval)
                        Term10 = Term10 + DubAmpl(j,m,b,e)*Retval
                     End Do
                  End Do
!!Having defined all the Terms in the amplitude equation, we combine them all together in the 'updated' excitation
!    amplitude:
                  UpdDubAmpl(i,j,a,b) = (ijab(i,j,a,b)+Term1-Term2-Term3&
&+Term4+Term5+Term6+Term7-Term8-Term9+Term10)/Dijab(i,j,a,b)
               End Do
            End Do
         End Do
      End Do
      AmplDiff = UpdDubAmpl-DubAmpl
      DubAmpl = UpdDubAmpl
!!!!!and thus the previous excitation amplitude is redefined as the updated value
!!Finally, the program recalculates the coupled-cluster Energy from the new amplitudes:
      Energy = 0
      Do b = 1, NV
         Do a = 1, NV
            Do j = 1, NO
               Do i = 1, NO
                  Energy = Energy + .25*DubAmpl(i,j,a,b)*ijab(i,j,a,b)
               End Do
            End Do
         End Do
      End Do
!!and checks for Convergence based on both the Energy dIfference (between current Energy and previous Energy), as well
!    as the dIfference in the root-mean squared (RMS) value of the excitation amplitudes:
      Call AmplConv(ConvA)
      ConvE = Energy - PrevEnergy
      ConvE = Sqrt(ConvE**2)
      If ((ConvA <= 1.0D-5) .and. (ConvE <= 1.0D-7)) then
         aaa = 1
      End If
!!!!!If Convergence has been achieved, the program signals the End of the iterative loop and prints the final Energy;
!         otherwise, it cycles back through the next iteration
      PrevEnergy = Energy
      Counter = Counter + 1
      Print*, "The energy of iteration ", Counter, "is", Energy
   End Do
   print*, "The final Energy is: ", Energy
End Subroutine CCD
!!The following blocks are each of the Subroutines/Terms incorporated into the overarching amplitude equation.
!
!!The F matrices correspond to their namesake and are defined by the nature of their indices - e.g. 'FV'
!         is the F matrix defined only through the virtual orbitals of the excitation.
Subroutine FV(b, e, Retval)
   Use Global
   Implicit Integer (f-o)
   Integer, intent(in) :: b, e
   Double Precision, intent(out) :: Retval
   Double Precision Temp, TaoT
   Integer NO, NV
   Common NO, NV
   Temp = 0
   Do f = 1, NV
      Do n = 1, NO
         Do m = 1, NO
            Temp = Temp + DubAmpl(m,n,b,f)*ijab(m,n,e,f)
         End Do
      End Do
   End Do
   Retval = Temp * (-.5)
End Subroutine FV

Subroutine FO(m, j, Retval)
   Use Global
   Implicit None
   Integer, intent(in) :: j, m
   Double Precision, intent(out) :: Retval
   Double Precision Temp, TaoT
   Integer n, e, f, NO, NV
   Common NO, NV
   Temp = 0
   Do f = 1, NV
      Do e = 1, NV
         Do n = 1, NO
            Temp = Temp + DubAmpl(j,n,e,f)*ijab(m,n,e,f)
         End Do
      End Do
   End Do
   Retval = Temp * .5
End Subroutine FO
!
!!The W matrices similarly correspond to their namesake and are defined by the nature of their indices - e.g. 'WV'
!         is the W matrix defined only through the virtual orbitals of the excitation.
Subroutine WO(m,n,i,j,Retval)
   Use Global
   Implicit None
   Integer, intent(in) :: m, n, i, j
   Double Precision, intent(out) :: Retval
   Double Precision Temp, Tao
   Integer e, f, NO, NV
   Common NO, NV
   Temp = 0
   Do f = 1, NV
      Do e = 1, NV
         Temp = Temp +.25*DubAmpl(i,j,e,f)*ijab(m,n,e,f)
      End Do
   End Do
   Retval = Temp + ijkl(m,n,i,j)
End Subroutine WO
!
Subroutine WV(a,b,e,f,Retval)
   Use Global
   Implicit None
   Integer, intent(in) :: a, b, e, f
   Double Precision, intent(out) :: Retval
   Double Precision Temp, Tao
   Integer m, n, NO, NV
   Common NO, NV
   Temp = 0
   Do n = 1, NO
      Do m = 1, NO
         Temp = Temp +.25*DubAmpl(m,n,a,b)*ijab(m,n,e,f)
      End Do
   End Do
   Retval = Temp + abcd(a,b,e,f)
End Subroutine WV
!
Subroutine WB(m,b,e,j,Retval)
   Use Global
   Implicit None
   Integer, intent(in) :: m, b, e, j
   Double Precision, intent(out) :: Retval
   Double Precision Temp, Tao
   Integer n, f, NO, NV
   Common NO, NV
   Temp = 0
   Do f = 1, NV
      Do n = 1, NO
         Tao = (-.5)*DubAmpl(j,n,f,b)! + sampl(j,f)*sampl(n,b)
         Temp = Temp + Tao*ijab(m,n,e,f)
      End Do
   End Do
   Retval = Temp + (-iajb(m,b,j,e))
End Subroutine WB
!!The 'AmplConv' Subroutine is a simple method of computing the RMS of the excitation Amplitudes in order to obtain a
!    value ('Retval', passed into the main program as 'ConvA') with which to check Convergence.
Subroutine AmplConv(Retval)
   Use Global
   Implicit None
   Double Precision, Intent(out) :: Retval
   Integer a, b, i, j, NO, NV
   Double Precision Temp, m
   Common NO, NV
   m = NO * NO * NV * NV
   Temp = 0
   Do b = 1, NV
      Do a = 1, NV
         Do j = 1, NO
            Do i = 1, NO
               Temp = Temp + AmplDiff(i,j,a,b)**2
            End Do
         End Do
      End Do
   End Do
   Temp = Temp / m
   Retval = sqrt(Temp)
End Subroutine AmplConv

max_issues_repo_name: caricato-ku/CHEM-Course-Code

max_forks_repo_name: caricato-ku/CHEM-Course-Code


text:
PROGRAM hello_world

  USE omp_lib

  PRINT *, 'Hello world from the main program'

  !$OMP PARALLEL

    PRINT '(A,I0)','Hello world! I am processor ', omp_get_thread_num()

  !$OMP END PARALLEL

END PROGRAM hello_world

max_stars_repo_name: WarwickRSE/OpenMPExamples

max_issues_repo_name: WarwickRSE/OpenMPExamples

max_forks_repo_name: WarwickRSE/OpenMPExamples


text:
 David Nguyens Biography:

Name: David Hoang Nguyen
Occupation: Freshman student at University of California, Davis
Phone: 408 239 6389
Resident: 492 La Rue Rd. 311A Davis, CA 95616
David Nguyen is a freshman student at University of California, Davis.  He typically enjoys watching movies, hanging out with friends, shooting pool, listen and play music, working out; he is up for anything that is fun.  David is currently a Computer Science major but intends to change it to something in the Biological Science field, perhaps something that relates to pharmaceuticals.
David Nguyen  Biography:
•	Born in San Jose, CA.
•	Graduated from Del Mar High School.

Plays the guitar/piano.

•	Enjoys watching movies.
•	Enjoys working out and playing sports.
•	Toured many regions of Asia.

20080515 18:31:05 nbsp welcome to the wiki David! Users/MattHh

20080516 10:07:48 nbsp Howdy, and welcome to the wiki!  Did you take that photo of La Rue Field yourself, or did you get it from a website?  If you dont have permission to reproduce it, it Copyrights cant be used on the wiki.  Great writeup, though! Users/JabberWokky


text:
!     This version solves the Schroedinger equation for the
!     one-dimensional
!     harmonic oscillator using matrix diagonalization

!     declaration of constants used in the calculations
!     define variables, where we have chosen to use atomic units for
!     m=c=hbar=k=1.

MODULE parameters
  USE constants
  !  the step
  REAL(DP), PUBLIC :: step
  !  r min
  REAL(DP), PUBLIC, PARAMETER :: rmin=-10.0_dp
  !  r max
  REAL(DP), PUBLIC, PARAMETER :: rmax=10.0_dp
  !  number of steps from rmin to rmin to rmax
  INTEGER, PUBLIC, PARAMETER ::  max_step=1000
  !  size of matrix to be diagonalized, do not need points at rmin and rmax
  INTEGER, PUBLIC, PARAMETER ::  n=998
END MODULE parameters

!  main program starts here

PROGRAM schroedinger_equation
  USE constants
  USE parameters
  USE f90library
  IMPLICIT NONE
  REAL(DP) :: potential, const1, const2
  REAL(DP), DIMENSION(:), ALLOCATABLE :: w, r
  REAL(DP), ALLOCATABLE, DIMENSION(:) :: e, d
  REAL(DP), ALLOCATABLE, DIMENSION(:,:) :: z
  INTEGER :: i

  !  reserve place in memory for the arrays w, r, e, d and z

  ALLOCATE ( e(n) , d(n))
  ALLOCATE ( w(0:max_step), r(0:max_step))

  !  define the step size

  step=(rmax-rmin)/FLOAT(max_step)

  !  define constants for the matrix to be diagonalized

  const1=2./(step*step)
  const2=-1./(step*step)

  !     set up r and the function w for energy =0
  !     w corresponds then to the potential
  !  values at

  DO i=0, max_step
     r(i) = rmin+i*step
     w(i) = potential(r(i))
  ENDDO

  !     setup the diagonal d and the non-diagonal part e  of
  !     the  tri-diagonal matrix to be diagonalized

  d(1:n)=const1+w(1:n)  ;  e(1:n)=const2

  !  allocate space for eigenvector info

  ALLOCATE ( z(n,n) )

  !  obtain the eigenvalues

  CALL tqli(d,e,n,z)

  !  sort eigenvalues as an ascending series

  CALL eigenvalue_sort(d,n)

  !  write out

  WRITE(6,*) ' Rmin and Rmax ', rmin, rmax
  WRITE(6,*) ' Number of steps -1 ', n
  WRITE(6,*) ' The 5 lowest eigenvalues '
  WRITE(6,'(5F12.6)') d(1:5)

  !  find the wave function for the lowest lying state

  CALL plot(d(1), w, r)

  !  deallocate space in memory

  DEALLOCATE (z)
  DEALLOCATE ( w, r, d, e)

END PROGRAM schroedinger_equation

!  Harmonic oscillator potential

REAL(DP) FUNCTION potential(x)
  USE constants
  IMPLICIT NONE
  REAL(DP), INTENT(IN) :: x

  potential=x*x

END FUNCTION  potential

!
!           Sort the eigenvalues as
!           ascending series
!

SUBROUTINE eigenvalue_sort(e,n)
  USE constants
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: n
  INTEGER :: i, j
  REAL(DP), DIMENSION(n), INTENT(INOUT) :: e
  REAL(DP) :: temp

  DO i = 1, n
     DO j =i, n
        IF( ABS(e(i)) >   ABS(e(j)) ) THEN
           ! change the energy
           temp=e(i)
           e(i)=e(j)
           e(j)=temp
        ENDIF
     ENDDO
  ENDDO

END SUBROUTINE eigenvalue_sort

!  This function plots the wave function at the final energy

SUBROUTINE plot(energy, w ,r )
  USE constants
  USE parameters
  IMPLICIT NONE
  INTEGER :: i, match
  REAL(DP), DIMENSION(n) :: w, r
  REAL(DP), DIMENSION(:), ALLOCATABLE :: wf
  REAL(DP) :: energy, wwf, fac, norm, pi

  ALLOCATE ( wf(0:max_step))
  pi=ACOS(-1.) ; pi=SQRT(pi)

  !  add the chosen energy to the potential

  w=(w-energy)*step*step+2

  !  Boundary conditions
  !  value of the wave function at max step

  wf(max_step)=0.

  !  here we use the analytical form for the wave function for max step -1
  !  (taylor expansion of EXP(-r^2/2))

  wf(max_step-1)=step**2/2.

  !  start from right and move towards rmin

  OUTER_TO_INNER : DO
     DO i=max_step-2, 1, -1
        wf(i)=w(i+1)*wf(i+1)-wf(i+2)
        !  find point were the outer part starts turning in order to determine matching point
        IF ( wf(i) <=  wf(i+1) ) EXIT OUTER_TO_INNER
     ENDDO
  ENDDO OUTER_TO_INNER
  match = i+1
  wwf=wf(match)

  !  start from rmin  and move towards the matching point

  wf(0)=0.
  wf(1)=step**2/2.
  DO i=2, match
     wf(i)=w(i-1)*wf(i-1)-wf(i-2)
  ENDDO

  !  scale the wave function below the matching point

  fac=wwf/wf(match)
  wf(0:match)=wf(0:match)*fac

  !  normalize the wave function

  norm=step*SUM(wf*wf)
  IF ( norm > 0. ) norm = 1./SQRT(norm)
  wf=wf*norm
  WRITE(6,*) norm
  !  write out of numerically calculated and analytical expression for the wave functions

  WRITE(6,'(E12.6,2X,E12.6,2X,E12.6)')( r(i), wf(i), EXP(-(r(i)**2)/2.0_dp)/SQRT(pi),i=1, max_step)
  DEALLOCATE ( wf)

END SUBROUTINE plot

max_stars_repo_name: kimrojas/ComputationalPhysicsMSU


text:
      subroutine sourceph(xp,yp,zp,nxp,nyp,nzp,
     +                    sint,cost,sinp,cosp,phi,fi,fq,fu,fv,
     +                    xmax,ymax,zmax,twopi,
     +                    xcell,ycell,zcell,nxg,nyg,nzg,iseed)

      implicit none

      include 'photon.txt'

      integer xcell,ycell,zcell,nxg,nyg,nzg,iseed
      real*8 xmax,ymax,zmax,twopi
      real ran2

c***** emit photon isotropically from origin
      xp=0.
      yp=0.
      zp=0.

      cost=2.*ran2(iseed)-1.
      sint=(1.-cost*cost)
      if(sint.le.0.)then
        sint=0.
      else
        sint=sqrt(sint)
      endif

      xp=0.
      yp=0.
      zp=0.

      phi=twopi*ran2(iseed)
      cosp=cos(phi)
      sinp=sin(phi)

c***** Set photon direction cosines for direction of travel *********
      nxp=sint*cosp
      nyp=sint*sinp
      nzp=cost

c***** Set Stokes fluxes ********************************************
      fi=1.
      fq=0.
      fu=0.
      fv=0.

c*************** Linear Grid *************************
      xcell=int(nxg*(xp+xmax)/(2.*xmax))+1
      ycell=int(nyg*(yp+ymax)/(2.*ymax))+1
      zcell=int(nzg*(zp+zmax)/(2.*zmax))+1
c*****************************************************

      return
      end


text:
MODULE MATRICES

  USE CONSTANTS

  USE GREEN_RANKINE
  USE GREEN_WAVE

  IMPLICIT NONE

CONTAINS

  SUBROUTINE ADD_RANKINE_PART_TO_THE_MATRICES(                        &
      nb_faces_1,                                                     &
      centers_1, normals_1,                                           &
      nb_vertices_2, nb_faces_2,                                      &
      vertices_2, faces_2, centers_2, normals_2, areas_2, radiuses_2, &
      coeff,                                                          &
      S, K)

    INTEGER,                                     INTENT(IN) :: nb_faces_1, nb_faces_2, nb_vertices_2
    REAL(KIND=PRE), DIMENSION(nb_faces_1, 3),    INTENT(IN) :: centers_1, normals_1
    REAL(KIND=PRE), DIMENSION(nb_vertices_2, 3), INTENT(IN) :: vertices_2
    INTEGER,        DIMENSION(nb_faces_2, 4),    INTENT(IN) :: faces_2
    REAL(KIND=PRE), DIMENSION(nb_faces_2, 3),    INTENT(IN) :: centers_2, normals_2
    REAL(KIND=PRE), DIMENSION(nb_faces_2),       INTENT(IN) :: areas_2, radiuses_2

    REAL(KIND=PRE), INTENT(IN) :: coeff

    COMPLEX(KIND=PRE), DIMENSION(nb_faces_1, nb_faces_2), INTENT(INOUT) :: S
    COMPLEX(KIND=PRE), DIMENSION(nb_faces_1, nb_faces_2), INTENT(INOUT) :: K

    ! Local variables
    INTEGER :: I, J
    REAL(KIND=PRE)               :: SP1
    REAL(KIND=PRE), DIMENSION(3) :: VSP1

    DO I = 1, nb_faces_1
      !$OMP PARALLEL DO PRIVATE(J, SP1, VSP1)
      DO J = 1, nb_faces_2

        CALL COMPUTE_INTEGRAL_OF_RANKINE_SOURCE( &
          centers_1(I, :),                       &
          vertices_2(faces_2(J, :), :),          &
          centers_2(J, :),                       &
          normals_2(J, :),                       &
          areas_2(J),                            &
          radiuses_2(J),                         &
          SP1, VSP1                              &
          )

        ! Store into influence matrix
        S(I, J) = S(I, J) - coeff * SP1/(4*PI)                                ! Green function
        K(I, J) = K(I, J) - coeff * DOT_PRODUCT(normals_1(I, :), VSP1)/(4*PI) ! Gradient of the Green function

      END DO
      !$OMP END PARALLEL DO
    END DO

  END SUBROUTINE

  ! =====================================================================

  SUBROUTINE ADD_WAVE_PART_TO_THE_MATRICES(  &
      nb_faces_1, centers_1, normals_1,  &
      nb_faces_2, nb_quad_points,        &
      quad_points, quad_weights,         &
      wavenumber, depth,                 &
      XR, XZ, APD,                       &
      NEXP, AMBDA, AR,                   &
      coeff,                             &
      same_body,                         &
      S, K)

    ! Mesh data
    INTEGER,                                  INTENT(IN) :: nb_faces_1, nb_faces_2
    REAL(KIND=PRE), DIMENSION(nb_faces_1, 3), INTENT(IN) :: normals_1, centers_1

    INTEGER,                                                  INTENT(IN) :: nb_quad_points
    REAL(KIND=PRE), DIMENSION(nb_faces_2, nb_quad_points, 3), INTENT(IN) :: quad_points
    REAL(KIND=PRE), DIMENSION(nb_faces_2, nb_quad_points),    INTENT(IN) :: quad_weights

    REAL(KIND=PRE),                           INTENT(IN) :: wavenumber, depth

    ! Tabulated integrals
    REAL(KIND=PRE), DIMENSION(328),           INTENT(IN) :: XR
    REAL(KIND=PRE), DIMENSION(46),            INTENT(IN) :: XZ
    REAL(KIND=PRE), DIMENSION(328, 46, 2, 2), INTENT(IN) :: APD

    ! Prony decomposition for finite depth
    INTEGER,                                  INTENT(IN) :: NEXP
    REAL(KIND=PRE), DIMENSION(NEXP),          INTENT(IN) :: AMBDA, AR

    REAL(KIND=PRE), INTENT(IN) :: coeff

    ! Trick to save some time
    LOGICAL,                                  INTENT(IN) :: same_body

    ! Output
    COMPLEX(KIND=PRE), DIMENSION(nb_faces_1, nb_faces_2), INTENT(INOUT) :: S
    COMPLEX(KIND=PRE), DIMENSION(nb_faces_1, nb_faces_2), INTENT(INOUT) :: K

    ! Local variables
    INTEGER                         :: I, J, Q
    COMPLEX(KIND=PRE)               :: SP2
    COMPLEX(KIND=PRE), DIMENSION(3) :: VSP2_SYM, VSP2_ANTISYM

    IF ((SAME_BODY) .AND. (nb_quad_points == 1)) THEN
      ! If we are computing the influence of some cells upon themselves, the resulting matrices have some symmetries.
      ! This is due to the symmetry of the Green function, and the way the integral on the face is approximated.
      ! (More precisely, the Green function is symmetric and its derivative is the sum of a symmetric part and an anti-symmetric
      ! part.)

      DO I = 1, nb_faces_1
        !$OMP PARALLEL DO PRIVATE(J, SP2, VSP2_SYM, VSP2_ANTISYM)
        DO J = I, nb_faces_2

          IF (depth == INFINITE_DEPTH) THEN
            CALL WAVE_PART_INFINITE_DEPTH &
              (wavenumber,                &
              centers_1(I, :),            &
              quad_points(J, 1, :),       & ! centers_2(J, :),
              XR, XZ, APD,                &
              SP2, VSP2_SYM               &
              )
            VSP2_ANTISYM(:) = ZERO
          ELSE
            CALL WAVE_PART_FINITE_DEPTH   &
              (wavenumber,                &
              centers_1(I, :),            &
              quad_points(J, 1, :),       & ! centers_2(J, :),
              depth,                      &
              XR, XZ, APD,                &
              NEXP, AMBDA, AR,            &
              SP2, VSP2_SYM, VSP2_ANTISYM &
              )
          END IF

          S(I, J) = S(I, J) - coeff/(4*PI) * SP2 * quad_weights(J, 1)
          K(I, J) = K(I, J) - coeff/(4*PI) * &
            DOT_PRODUCT(normals_1(I, :), VSP2_SYM + VSP2_ANTISYM) * quad_weights(J, 1)

          IF (.NOT. I==J) THEN
            VSP2_SYM(1:2) = -VSP2_SYM(1:2)
            S(J, I) = S(J, I) - coeff/(4*PI) * SP2 * quad_weights(I, 1)
            K(J, I) = K(J, I) - coeff/(4*PI) * &
              DOT_PRODUCT(normals_1(J, :), VSP2_SYM - VSP2_ANTISYM) * quad_weights(I, 1)
          END IF

        END DO
        !$OMP END PARALLEL DO
      END DO

    ELSE
      ! General case: if we are computing the influence of a some cells on other cells, we have to compute all the coefficients.

      DO I = 1, nb_faces_1
        !$OMP PARALLEL DO PRIVATE(J, SP2, VSP2_SYM, VSP2_ANTISYM)
        DO J = 1, nb_faces_2
          DO Q = 1, nb_quad_points
            IF (depth == INFINITE_DEPTH) THEN
              CALL WAVE_PART_INFINITE_DEPTH &
                (wavenumber,                &
                centers_1(I, :),            &
                quad_points(J, Q, :),       &
                XR, XZ, APD,                &
                SP2, VSP2_SYM               &
                )
              VSP2_ANTISYM(:) = ZERO
            ELSE
              CALL WAVE_PART_FINITE_DEPTH   &
                (wavenumber,                &
                centers_1(I, :),            &
                quad_points(J, Q, :),       &
                depth,                      &
                XR, XZ, APD,                &
                NEXP, AMBDA, AR,            &
                SP2, VSP2_SYM, VSP2_ANTISYM &
                )
            END IF

            S(I, J) = S(I, J) - coeff/(4*PI) * SP2 * quad_weights(J, Q)
            K(I, J) = K(I, J) - coeff/(4*PI) * &
              DOT_PRODUCT(normals_1(I, :), VSP2_SYM + VSP2_ANTISYM) * quad_weights(J, Q)

          END DO
        END DO
        !$OMP END PARALLEL DO
      END DO
    END IF

  END SUBROUTINE

  ! =====================================================================

  SUBROUTINE BUILD_MATRICES(                          &
      nb_faces_1, centers_1, normals_1,               &
      nb_vertices_2, nb_faces_2, vertices_2, faces_2, &
      centers_2, normals_2, areas_2, radiuses_2,      &
      nb_quad_points, quad_points, quad_weights,      &
      wavenumber, depth,                              &
      coeffs,                                         &
      XR, XZ, APD,                                    &
      NEXP, AMBDA, AR,                                &
      same_body,                                      &
      S, K)

    ! Mesh data
    INTEGER,                                     INTENT(IN) :: nb_faces_1, nb_faces_2, nb_vertices_2
    REAL(KIND=PRE), DIMENSION(nb_faces_1, 3),    INTENT(IN) :: centers_1, normals_1
    REAL(KIND=PRE), DIMENSION(nb_vertices_2, 3), INTENT(IN) :: vertices_2
    INTEGER,        DIMENSION(nb_faces_2, 4),    INTENT(IN) :: faces_2
    REAL(KIND=PRE), DIMENSION(nb_faces_2, 3),    INTENT(IN) :: centers_2, normals_2
    REAL(KIND=PRE), DIMENSION(nb_faces_2),       INTENT(IN) :: areas_2, radiuses_2

    INTEGER,                                                  INTENT(IN) :: nb_quad_points
    REAL(KIND=PRE), DIMENSION(nb_faces_2, nb_quad_points, 3), INTENT(IN) :: quad_points
    REAL(KIND=PRE), DIMENSION(nb_faces_2, nb_quad_points),    INTENT(IN) :: quad_weights

    LOGICAL,                                  INTENT(IN) :: same_body

    REAL(KIND=PRE),                           INTENT(IN) :: wavenumber, depth

    REAL(KIND=PRE), DIMENSION(3) :: coeffs

    ! Tabulated integrals
    REAL(KIND=PRE), DIMENSION(328),           INTENT(IN) :: XR
    REAL(KIND=PRE), DIMENSION(46),            INTENT(IN) :: XZ
    REAL(KIND=PRE), DIMENSION(328, 46, 2, 2), INTENT(IN) :: APD

    ! Prony decomposition for finite depth
    INTEGER,                                  INTENT(IN) :: NEXP
    REAL(KIND=PRE), DIMENSION(NEXP),          INTENT(IN) :: AMBDA, AR

    ! Output
    COMPLEX(KIND=PRE), DIMENSION(nb_faces_1, nb_faces_2), INTENT(OUT) :: S
    COMPLEX(KIND=PRE), DIMENSION(nb_faces_1, nb_faces_2), INTENT(OUT) :: K

    ! Local variables
    INTEGER :: I
    REAL(KIND=PRE), DIMENSION(nb_faces_1, 3) :: reflected_centers_1, reflected_normals_1

    !!!!!!!!!!!!!!!!!!!!
    !  Initialization  !
    !!!!!!!!!!!!!!!!!!!!

    S(:, :) = CMPLX(0.0, 0.0, KIND=PRE)
    K(:, :) = CMPLX(0.0, 0.0, KIND=PRE)

    !!!!!!!!!!!!!!!!!!
    !  Rankine part  !
    !!!!!!!!!!!!!!!!!!

    IF (coeffs(1) .NE. ZERO) THEN
      CALL ADD_RANKINE_PART_TO_THE_MATRICES(                            &
        nb_faces_1,                                                     &
        centers_1, normals_1,                                           &
        nb_vertices_2, nb_faces_2,                                      &
        vertices_2, faces_2, centers_2, normals_2, areas_2, radiuses_2, &
        coeffs(1),                                                      &
        S, K)
    END IF

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !  Reflected Rankine part  !
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!

    IF (coeffs(2) .NE. ZERO) THEN

      IF (depth == INFINITE_DEPTH) THEN
        ! Reflection through free surface
        reflected_centers_1(:, 1:2) = centers_1(:, 1:2)
        reflected_centers_1(:, 3)   = -centers_1(:, 3)
      ELSE
        ! Reflection through sea bottom
        reflected_centers_1(:, 1:2) = centers_1(:, 1:2)
        reflected_centers_1(:, 3)   = -centers_1(:, 3) - 2*depth
      END IF

      reflected_normals_1(:, 1:2) = normals_1(:, 1:2)
      reflected_normals_1(:, 3)   = -normals_1(:, 3)

      CALL ADD_RANKINE_PART_TO_THE_MATRICES(                            &
        nb_faces_1,                                                     &
        reflected_centers_1, reflected_normals_1,                       &
        nb_vertices_2, nb_faces_2,                                      &
        vertices_2, faces_2, centers_2, normals_2, areas_2, radiuses_2, &
        coeffs(2),                                                      &
        S, K)

    END IF

    !!!!!!!!!!!!!!!
    !  Wave part  !
    !!!!!!!!!!!!!!!

    IF (coeffs(3) .NE. ZERO) THEN

      CALL ADD_WAVE_PART_TO_THE_MATRICES(  &
        nb_faces_1, centers_1, normals_1,  &
        nb_faces_2, nb_quad_points,        &
        quad_points, quad_weights,         &
        wavenumber, depth,                 &
        XR, XZ, APD,                       &
        NEXP, AMBDA, AR,                   &
        coeffs(3),                         &
        same_body,                         &
        S, K)

    END IF

    !!!!!!!!!!!!!

    IF (SAME_BODY) THEN
      DO I = 1, nb_faces_1
        K(I, I) = K(I, I) + 0.5
      END DO
    END IF

    RETURN

  END SUBROUTINE

  ! =====================================================================

END MODULE MATRICES

max_stars_repo_name: mancellin/libDelhommeau

max_issues_repo_name: mancellin/libDelhommeau

max_forks_repo_name: mancellin/libDelhommeau


text:
module ExampleOptionalParameter
  ! use any module needed for the sort function(s)
  ! and all the interfaces needed to make the code work
  implicit none
contains

  subroutine sort_table(table, ordering, column, reverse)
    type(table_type), intent(inout) :: table
    integer, optional :: column
    logical, optional :: reverse
    optional :: ordering
    interface
       integer function ordering(a, b)
         type(table_element), intent(in) :: a, b
       end function ordering
    end interface

    integer :: the_column, i
    logical :: reversing
    type(table_row) :: rowA, rowB

    if ( present(column) ) then
       if ( column > get_num_of_columns(table) ) then
          ! raise an error?
       else
          the_column = column
       end if
    else
       the_column = 1   ! a default value, de facto
    end if

    reversing = .false.  ! default value
    if ( present(reverse) ) reversing = reverse

    do
       ! loops over the rows to sort... at some point, we need
       ! comparing an element (cell) of the row, with the element
       ! in another row; ... let us suppose rowA and rowB are
       ! the two rows we are considering
       ea = get_element(rowA, the_column)
       eb = get_element(rowB, the_column)
       if ( present(ordering) ) then
          if ( .not. reversing ) then
             if ( ordering(ea, eb) > 0 ) then
                ! swap the rowA with the rowB
             end if
          else   ! < instead of >
             if ( ordering(ea, eb) < 0 ) then
                ! swap the rowA with the rowB
             end if
          end if
       else
          if ( .not. reversing ) then
             if ( lexinternal(ea, eb) > 0 ) then
                ! swap the rowA with the rowB
             end if
          else   ! < instead of >
             if ( lexinternal(ea, eb) < 0 ) then
                ! swap the rowA with the rowB
             end if
          end if
       end if
       ! ... more of the sorting algo ...
       ! ... and rows traversing ... (and an exit condition of course!)
    end do

  end subroutine sort_table

end module ExampleOptionalParameter

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


! This subroutine is used to caculate Hamiltonian for
! slab system .

! History
!        4/18/2010 by Quansheng Wu
! Copyright (c) 2010 QuanSheng Wu. All rights reserved.
  subroutine ham_slab(k,Hamk_slab)

     use para
     implicit none

     ! loop index
     integer :: i1, i2, i, j
     integer :: ir

     ! wave vector in 2d
     real(Dp), intent(inout) :: k(2)

     ! Hamiltonian of slab system
     complex(Dp),intent(out) ::Hamk_slab(Num_wann*nslab,Num_wann*nslab)

     ! the factor 2 is induced by spin
     complex(Dp), allocatable :: Hij(:, :, :)

     allocate( Hij(-ijmax:ijmax,Num_wann,Num_wann))
     call ham_qlayer2qlayer2(k,Hij)

     Hamk_slab=0.0d0
     ! i1 column index
     do i1=1, nslab
        ! i2 row index
        do i2=1, nslab
          if (abs(i2-i1).le.ijmax)then
            Hamk_slab((i2-1)*Num_wann+1:(i2-1)*Num_wann+Num_wann,&
                      (i1-1)*Num_wann+1:(i1-1)*Num_wann+Num_wann )&
            = Hij(i2-i1,1:Num_wann,1:Num_wann)
          endif
        enddo ! i2
     enddo ! i1

     ! check hermitcity

     do i1=1,nslab*Num_wann
     do i2=1,nslab*Num_wann
        if(abs(Hamk_slab(i1,i2)-conjg(Hamk_slab(i2,i1))).ge.1e-6)then
          write(stdout,*)'there are something wrong with Hamk_slab'
          stop
        endif
     enddo
     enddo

  return
  end subroutine ham_slab

  ! This subroutine is used to caculate Hamiltonian for
  ! slab system .
  !> for slab with in-plane magnetic field
  !> the magnetic phase are chosen like this
  !> phi_ij= alpha*[By*(xj-xi)*(zi+zj)-Bx*(yj-yi)*(zi+zj)]
  !> x, z are in unit of Angstrom, Bx, By are in unit of Tesla
  !> History :
  !        9/21/2015 by Quansheng Wu @ETH Zurich
  subroutine ham_slab_parallel_B(k,Hamk_slab)

     use para
     implicit none

     ! loop index
     integer :: i1, i2, i, j

     ! wave vector in 2d
     real(Dp), intent(inout) :: k(2)

     ! loop index
     integer :: iR

     ! index used to sign irvec
     integer :: ia,ib,ic
     integer :: ia1, ia2

     integer :: istart1, istart2
     integer :: iend1, iend2

     integer :: inew_ic

     !> nwann= Num_wann/2
     integer :: nwann

     integer, allocatable :: orbital_start(:)

     ! new index used to sign irvec
     real(dp) :: new_ia,new_ib,new_ic

     ! wave vector k times lattice vector R
     real(Dp) :: kdotr
     real(dp) :: phase
     complex(dp) :: ratio
     complex(dp) :: fac

     real(dp) :: R(3)
     real(dp) :: Rp1(3)
     real(dp) :: Rp2(3)
     real(dp) :: R1(3)
     real(dp) :: R2(3)
     real(dp) :: Ri(3)
     real(dp) :: Rj(3)
     real(dp) :: tau1(3)
     real(dp) :: tau2(3)

     ! Hamiltonian of slab system
     complex(Dp),intent(out) ::Hamk_slab(Num_wann*nslab,Num_wann*nslab)

     nwann= Num_wann/2
     allocate( orbital_start(Num_atoms+ 1))
     orbital_start= 0
     orbital_start(1)= 1
     do ia=1, Num_atoms
        orbital_start(ia+1)= orbital_start(ia)+ nprojs(ia)
     enddo

     Hamk_slab=0.0d0
     ! i1 column index
     do i1=1, Nslab
        ! i2 row index
        do i2=1, Nslab
           !> sum over R points to get H(k1, k2)
           do iR=1, Nrpts
              ia=irvec(1,iR)
              ib=irvec(2,iR)
              ic=irvec(3,iR)

              !> new lattice
              call latticetransform(ia, ib, ic, new_ia, new_ib, new_ic)

              inew_ic= int(new_ic)
              if (inew_ic /= (i2-i1)) cycle

              !> exp(i k.R)
              kdotr= k(1)*new_ia+ k(2)*new_ib
              ratio= cos(2d0*pi*kdotr)+ zi*sin(2d0*pi*kdotr)

              R1= (i1-1)*Ruc_new
              R2= new_ia*Rua_new+ new_ib*Rub_new+ (i2-1)*Ruc_new

              do ia1=1, Num_atoms
              do ia2=1, Num_atoms
                 R1= Atom_position(:, ia1)
                 R2= Atom_position(:, ia2)
                 call rotate(R1, tau1)
                 call rotate(R2, tau2)

                 Ri= Rp1+ tau1
                 Rj= Rp2+ tau2

                 phase= alpha*By*(Rj(3)+Ri(3))*(Rj(1)-Ri(1))  &
                      - alpha*Bx*(Rj(3)+Ri(3))*(Rj(2)-Ri(2))
                 fac= cos(phase)+ zi*sin(phase)

                !write(*, '(a, 4i5   )') 'iR, ia ib ic', ir, ia, ib, ic
                !write(*, '(a, 4f10.5)') '            ', new_ia, new_ib, new_ic
                !write(*, '(a, 3f10.5)') 'Ri', Ri
                !write(*, '(a, 3f10.5)') 'Rj', Rj
                !write(*, '(a, 3f10.5)') 'phase', phase

                 istart1= (i2-1)*Num_wann+ orbital_start(ia1)
                 iend1= (i2-1)*Num_wann+ orbital_start(ia1+1)- 1
                 istart2= (i1-1)*Num_wann+ orbital_start(ia2)
                 iend2= (i1-1)*Num_wann+ orbital_start(ia2+1)- 1

                 Hamk_slab( istart1:iend1, istart2:iend2) &
                 = Hamk_slab( istart1:iend1, istart2:iend2) &
                 + HmnR( istart1- (i2-1)*Num_wann:iend1- (i2-1)*Num_wann, &
                 istart2- (i1-1)*Num_wann:iend2- (i1-1)*Num_wann, iR)*ratio/ndegen(iR)* fac

                 !> there is soc term in the hr file
                 if (soc>0) then
                    istart1= (i2-1)*Num_wann+ orbital_start(ia1) + Nwann
                    iend1= (i2-1)*Num_wann+ orbital_start(ia1+1)- 1 + Nwann
                    istart2= (i1-1)*Num_wann+ orbital_start(ia2)
                    iend2= (i1-1)*Num_wann+ orbital_start(ia2+1)- 1

                    Hamk_slab( istart1:iend1, istart2:iend2) &
                    = Hamk_slab( istart1:iend1, istart2:iend2) &
                    + HmnR( istart1- (i2-1)*Num_wann:iend1- (i2-1)*Num_wann, &
                    istart2- (i1-1)*Num_wann:iend2- (i1-1)*Num_wann, iR)*ratio/ndegen(iR)* fac

                    istart1= (i2-1)*Num_wann+ orbital_start(ia1)
                    iend1= (i2-1)*Num_wann+ orbital_start(ia1+1)- 1
                    istart2= (i1-1)*Num_wann+ orbital_start(ia2) + Nwann
                    iend2= (i1-1)*Num_wann+ orbital_start(ia2+1)- 1 + Nwann

                    Hamk_slab( istart1:iend1, istart2:iend2) &
                    = Hamk_slab( istart1:iend1, istart2:iend2) &
                    + HmnR( istart1- (i2-1)*Num_wann:iend1- (i2-1)*Num_wann, &
                    istart2- (i1-1)*Num_wann:iend2- (i1-1)*Num_wann, iR)*ratio/ndegen(iR)* fac

                    istart1= (i2-1)*Num_wann+ orbital_start(ia1) + Nwann
                    iend1= (i2-1)*Num_wann+ orbital_start(ia1+1)- 1 + Nwann
                    istart2= (i1-1)*Num_wann+ orbital_start(ia2) + Nwann
                    iend2= (i1-1)*Num_wann+ orbital_start(ia2+1)- 1 + Nwann

                    Hamk_slab( istart1:iend1, istart2:iend2) &
                    = Hamk_slab( istart1:iend1, istart2:iend2) &
                    + HmnR( istart1- (i2-1)*Num_wann:iend1- (i2-1)*Num_wann, &
                    istart2- (i1-1)*Num_wann:iend2- (i1-1)*Num_wann, iR)*ratio/ndegen(iR)* fac
                 endif ! soc
              enddo ! ia2
              enddo ! ia1
           enddo ! iR
        enddo ! i2
     enddo ! i1

     !> check hermitcity
     do i1=1,nslab*Num_wann
     do i2=1,nslab*Num_wann
        if(abs(Hamk_slab(i1,i2)-conjg(Hamk_slab(i2,i1))).ge.1e-6)then
          write(stdout,*)'there are something wrong with Hamk_slab'
          stop
        endif
     enddo
     enddo

  return
  end subroutine ham_slab_parallel_B


text:
2205 Haring, located at the far west end of Haring Hall is probably the least comfortable Campus Lecture Halls lecture hall on campus. Built when students were significantly smaller than they are now, the seats in Haring are far too small for the average student. Wise students will arrive early and get seats in the second to front row, granting themselves some leg room.

Amazingly, I always tend to fall asleep in there.  Users/RohiniJasavala Ro

Summer Advising 2004 was held in there.  It was so hard to complete the chemistry and precalculus diagnostic tests in those horrible seats.  Users/PatrickSing

If you dont like the smallness of the seats here, go find a seat at the ARC so you can stop being obese. Nyah!  Users/Yesitsme


text:
Subroutine wida1(dmass, rhomp, wa1, iseed)
  Save
  pimass = 0.137265
  coupa = 14.8
  rhomax = dmass - pimass - 0.02
  If (rhomax<=0) Then
     rhomp = 0.
     wa1 = -10.
  End If
  icount = 0
711 rhomp = rhomas(rhomax, iseed)
  icount = icount + 1
  If (dmass<=(pimass+rhomp)) Then
     If (icount<=100) Then
        Goto 711
     Else
        rhomp = 0.
        wa1 = -10.
        Return
     End If
  End If
  qqp2 = (dmass**2-(rhomp+pimass)**2)*(dmass**2-(rhomp-pimass)**2)
  qqp = sqrt(qqp2)/(2.0*dmass)
  epi = sqrt(pimass**2+qqp**2)
  erho = sqrt(rhomp**2+qqp**2)
  epirho = 2.0*(epi*erho+qqp**2)**2 + rhomp**2*epi**2
  wa1 = coupa**2*qqp*epirho/(24.0*3.1416*dmass**2)
  Return
End Subroutine wida1

max_stars_repo_name: xiaohaijin/AMPT

max_issues_repo_name: xiaohaijin/AMPT

max_forks_repo_name: xiaohaijin/AMPT


text: Jeremy Day is one of several Edward Jones representatives in town.


text:
      subroutine best_xyz(bestx,besty,bestz,cal,srtd,hes,Tend,
     &     vtx,montecarlo)
C
C best_xyz returns the best reconstructed position for the scattered
c     electron given the CAL, SRTD, HES and track endpoint for the
c     electron and the measured vertex.
C
C  For Rcal electrons (Z<-140) the Z is projected onto the CAL plane.
C
C
C INPUT:
C    Real cal(4)
C         cal(1) = x cal from elecpo
C         cal(2) = y cal from elecpo
C         cal(3) = z cal from elecpo ( not used)
C         cal(4) = electron energy (not used by best_theta)
C
C    Real srtd(4)
C         srtd(1) = x from srtdelec
C         srtd(2) = y from srtdelec
C         srtd(3) = z from srtdelec (not used)
C         srtd(4) = srtd mips
C                   (if mips is in the range 0.4 to 1000 then srtd can be used)
C                   (NO srtd mips are returned by srtdelec if less than 0.4)
C
C    Real hes(4)
C         hes(1) = x from hesclu
C         hes(2) = y from hesclu
C         hes(3) = z from hesclu (not used)
C         hes(4) = hes mips
C                   (if mips is in the range 5 to 1000 then hes can be used)
C                   (NO HES mips are returned by hesclu if less than 5 mips)
C    Real Tend(5)
C         tend(1) = x from Track
C         tend(2) = y from Track
C         tend(3) = z from track (not used)
C         tend(4) = Track Momentum
C                   (if Momentum > 5 then Tend can be used)
C         tend(5) = track DCA
C                   (DCA must be < 10 for this to be used)
C
C    Real vtx(3) = vertex (x,y,z) position
C
C    Logical Montecarlo = True if montecarlo
C
C
C Author: Mike Wodarczyk  13-June-1998   Version 1.0
C
C
      Implicit None

      Real cal(4),srtd(4),hes(4),tend(5),vtx(3)
      Logical Montecarlo
      Real bestx, besty,bestz

      real RCAL_FACE
      parameter ( RCAL_FACE   =   -152.13 )
      real BCAL_FACE
      parameter ( BCAL_FACE   =   125.669 )
      real FCAL_FACE
      parameter ( FCAL_FACE   =   226.13  )

      real dataRCALz,mcRCALz, rcalz
      parameter ( dataRCALz =   -153.88 )
      parameter ( mcRCALz   =   -152.13 )
      real datahesz,mchesz, hesz
      parameter ( datahesz  =   -155.21 ) ! best position for no bias
      parameter ( mchesz    =   -154.14 ) ! mc hes Z
      real dataFhesz,mcFhesz, fhesz
      parameter ( dataFhesz  =   227.2 )
      parameter ( mcFhesz    =   227.2 )
      real datasrtdz, mcsrtdz, srtdz
      parameter ( datasrtdz =  -149.31 )  ! same as HES Z shift
      parameter ( mcsrtdz   =  -148.25 )  ! midpoint of 2 srtd planes

      integer whatToUse
      integer useCAL, useHES, useSRTD, useTrack
      parameter ( useCAL    = 1 )
      parameter ( useHES    = 2 )
      parameter ( useSRTD   = 3 )
      parameter ( useTrack = 4 )

      real calradius, radius
      real modpos
      real hesgap
C     leave a gap of 2cm on BOTH SIDES of Module crack
      parameter ( hesgap = 2.0 )
      real SRTDedge
C     Use SRTD away from the edge by SRTDedge cm
      parameter ( SRTDedge = 2.0 )

      real MAXMIP
      PARAMETER ( MAXMIP = 1000. ) ! stay away from high default values.

      if (montecarlo) then
         rcalz = mcrcalz
         hesz  = mchesz
         fhesz = mcfhesz
         srtdz  = mcsrtdz
      else
         rcalz = datarcalz
         hesz  = datahesz
         fhesz = datafhesz
         srtdz  = datasrtdz
      endif

      WhatToUse = useCAL  ! start by using cal
      calradius = sqrt(cal(1)**2+cal(2)**2)

C -------------------------------
C  See if we can use HES
C -------------------------------

C
C  IN RCAL use HES over CAL when not in module cracks and
C  the position is more than 20 cm from RCAL beampipe
C  (near the beampipe the hes diodes end and the SRTD
C  takes over anyway.
C
      if ( cal(3).lt.-140 .and. calRadius.gt.20.) then
C     in the RCAL more than 20 cm from beampipe.

         modpos = abs(cal(1)) - 20.33/2. + hesgap
         if ( modpos - int(modpos/20.33)*20.33 .gt. 2*hesgap .or.
     &        modpos.lt.0 ) then
C
C     The cal says that we are NOT IN THE MODULE CRACK so I expect there
C     to be a HES MIP.
C
            if (hes(4).ge.5..and.hes(4).lt.MAXMIP) then
               WhatToUse = UseHES
            endif
         endif ! not in crack
      endif ! in rcal
C
C  In FCAL use the hes whenever available
C
      if ( cal(3).gt.220.and.hes(4).gt.5.and.hes(4).lt.MAXMIP) then
         WhatToUse = UseHES
      endif

C -------------------------------
C  See if we can use SRTD
C
C  require an cal position to be in RCAL and the SRTD mips to be o.k.
C  then require that the electron is within SRTDedge cm of the srtd edge.
C
C -------------------------------
      if ( cal(3).lt.-140 .and.srtd(4).gt.0.4.and.srtd(4).lt.MAXMIP )
     $     then
         if (srtd(1).gt.(-34+srtdedge).and.srtd(1).lt.(34-srtdedge).and
     $        .(srtd(1).gt.(-10-srtdedge).and.srtd(2).gt.(-28+srtdedge)
     $        .or.srtd(1).le.(-10-srtdedge).and.srtd(2).gt.-40+srtdedge)
     $        .and.(srtd(1).ge.(10+srtdedge).and.srtd(2).lt.(40-srtdedge
     $        ).or.srtd(1).lt.(10+srtdedge).and.srtd(2).lt.28-srtdedge)
     $        )then
            WhatToUse = UseSRTD
         endif
      endif

C -------------------------------
C  See if we can use Track Endpoint
C
C  We must be more than 80 cm away from 0,0 for there to be
C  enough stereo hits for the track to give good resolution.
C -------------------------------
      if ( calRadius.gt.80.and.tend(4).gt.5.and.tend(5).lt.10) then
         WhatToUse = UseTrack
      endif

C
C  Start with Elecpo
C
      if ( WhatToUse .eq. UseCal ) then
         if ( cal(3).lt.-140 ) then
C     ...  USE SHIFTED RCAL
            bestx = cal(1)
            besty = cal(2)
            bestz = rcalz
         else
            bestx = cal(1)
            besty = cal(2)
            bestz = cal(3)
         endif
      endif

C
C  Calculate HES theta
C
      if ( WhatToUse .eq. UseHES ) then
         if ( cal(3).lt.-140 ) then
C     ...  USE SHIFTED RCAL
            bestx = hes(1)
            besty = hes(2)
            bestz = hesz
         else
            bestx = hes(1)
            besty = hes(2)
            bestz = fhesz
         endif
      endif

C
C  Calculate SRTD Theta
C
      if ( WhatToUse .eq. UseSRTD ) then
         bestx = srtd(1)
         besty = srtd(2)
         bestz = srtdz
      endif

C
C  Calculate VCTRHL Theta
C
      if ( WhatToUse .eq. UseTrack) then
         bestx = Tend(1)
         besty = Tend(2)
         bestz = Tend(3)
      endif

C
C  Now project all x,y,z onto planes at the RCAL FACE or Bcal cylinder
c     or FCAL face
C

      if ( bestz .lt. -140 ) then
C RCAL
         bestx = bestx * (RCAL_FACE - vtx(3)) / (bestz - vtx(3))
         besty = besty * (RCAL_FACE - vtx(3)) / (bestz - vtx(3))
         bestz = RCAL_FACE
      else if (bestz.lt.220 ) then
C BCAL
         radius = sqrt(bestx**2+besty**2)
         bestx = bestx*BCAL_FACE / radius
         besty = besty*BCAL_FACE / radius
         bestz = vtx(3) + (bestz-vtx(3))*BCAL_FACE / radius
      else
C FCAL
         bestx = bestx * (FCAL_FACE - vtx(3)) / (bestz - vtx(3))
         besty = besty * (FCAL_FACE - vtx(3)) / (bestz - vtx(3))
         bestz = FCAL_FACE
      endif

      return
      end

ext: fpp


text:
      SUBROUTINE eval_saddle_coils
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE stel_constants
      USE boundary, ONLY: nfp
      USE saddle_coils
      USE saddle_surface
      USE Vcoilpts
      USE Vwire
      USE coils
      IMPLICIT NONE
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: i, j, k, n, m, nc, np, ierr
      REAL(rprec), ALLOCATABLE, DIMENSION(:) :: uc, us
      REAL(rprec), ALLOCATABLE, DIMENSION(:) :: vc, vs
      REAL(rprec) :: u,  v,  s, ds, cmax, dlsq, csum, crad
      REAL(rprec), DIMENSION(2,3,5) :: x,  y,  z
      REAL(rprec), DIMENSION(2) :: c
      REAL(rprec) :: r, alf, bet
      LOGICAL :: la, lb
!-----------------------------------------------
      ALLOCATE (uc(0:nsad_u), us(0:nsad_u+4),
     1          vc(0:nsad_v), vs(0:nsad_v+4))

      alf = 1
      bet = 0
      la = .true.
      lb = .false.

      nc = nsad_coils_per_period
      ds = 1.0_dp/nwire
!     loop over unique coils
      DO n=1, nsmid
         DO k=0, nsad_u
            uc(k) = saddle(n)%u_c(k)
            us(k) = saddle(n)%u_s(k)
         END DO
         DO k=0, nsad_v
            vc(k) = saddle(n)%v_c(k)
            vs(k) = saddle(n)%v_s(k)
         END DO
         m = nc + 1 - n
!        loop over number of poloidal points
         cmax = 0
         csum = 0
         DO i=1, nwire+1
            s = (i-1)*ds
            CALL coil_curve( nsad_u, uc, us, nsad_v, vc, vs, la, lb,
     1         lctrlpt, lspline, nfp, numsurf_sad, m_sad, n_sad,
     2         rmn_sad, zmn_sad, nfils, s, alf, bet, deln, delt,
     3         u, v, r, c, x, y, z, ierr )
            cmax = MAX (cmax, c(1))
            csum = csum + c(1)
!           set coil x,y,z points for ALL field periods
            DO j=1, nfp
               np = (j-1)*nc
               u_sad(i,n + np) = u
               v_sad(i,n + np) = v + 1.0*(j-1)
!              stellarator symmetry
               u_sad(nwire + 2 - i,m + np) = 1.0 - u
               v_sad(nwire + 2 - i,m + np) = -v + 1.0*(j-1)
               DO k=1, nfils
                  x_sad(i,n + np,k) = x(1,j,k)
                  y_sad(i,n + np,k) = y(1,j,k)
                  z_sad(i,n + np,k) = z(1,j,k)
!                 stellarator symmetry
                  x_sad(i,m + np,k) = x(2,j,k)
                  y_sad(i,m + np,k) = y(2,j,k)
                  z_sad(i,m + np,k) = z(2,j,k)
               END DO
            END DO
         END DO
         rs_min(n) = 10000.0_dp
         IF(cmax .gt. 0.0_dp) rs_min(n) = 1.0_dp/cmax
         cs_sum(n) = csum*ds
      END DO

!     coil currents
      DO n=1, nsmid
         m = nc + 1 - n
         DO j=1, nfp
            np = (j-1)*nc
            c_sad(n + np) =  saddle(n)%current
            c_sad(m + np) = -saddle(n)%current
         END DO
      END DO

!     coil length
      DO n=1, nsmid
         ymin_sad(n) = 1000
         rmax_sad(n) = 0
         sad_length(n) = 0
         DO i=1, nwire
            dlsq = (x_sad(i+1,n,1) - x_sad(i,n,1))**2
     1           + (y_sad(i+1,n,1) - y_sad(i,n,1))**2
     2           + (z_sad(i+1,n,1) - z_sad(i,n,1))**2
            sad_length(n) = sad_length(n) + SQRT(dlsq)
            IF (ABS(y_sad(i,n,1)) .lt. ymin_sad(n))
     1         ymin_sad(n) = ABS(y_sad(i,n,1))
            crad = SQRT(x_sad(i,n,1)**2+y_sad(i,n,1)**2)
            IF (crad .gt. rmax_sad(n)) rmax_sad(n) = crad
         END DO
      END DO

      DEALLOCATE (uc, us, vc, vs)

      END SUBROUTINE eval_saddle_coils

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: joseluisvelasco/STELLOPT


text:
      INTEGER FUNCTION IDAMAX(N,X,INCX)
c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
c ALL RIGHTS RESERVED.
c Based on Government Sponsored Research NAS7-03001.
C>> 1994-11-11 IDAMAX  Krogh   Declared all vars.
C>> 1994-10-20 IDAMAX Krogh  Changes to use M77CON
C>> 1994-04-19 IDAMAX Krogh   Conveted to use generic intrinsics.
C>> 1987-12-09 IDAMAX Lawson  Initial code.
c--D replaces "?": I?AMAX
C
C     FIND SMALLEST INDEX OF MAXIMUM MAGNITUDE OF X.
C     IDAMAX =  FIRST I, I = 1 TO N, TO MINIMIZE  ABS(X(1-INCX+I*INCX))
C
      INTEGER N, INCX, NS, II, I
      DOUBLE PRECISION X(*),XMAX,XMAG
      IDAMAX = 0
      IF(N.LE.0) RETURN
      IDAMAX = 1
      IF(N.LE.1)RETURN
      IF(INCX.EQ.1)GOTO 20
C
C        CODE FOR INCREMENTS NOT EQUAL TO 1.
C
      XMAX = ABS(X(1))
      NS = N*INCX
      II = 1
          DO 10 I = 1,NS,INCX
          XMAG = ABS(X(I))
          IF(XMAG.LE.XMAX) GO TO 5
          IDAMAX = II
          XMAX = XMAG
    5     II = II + 1
   10     CONTINUE
      RETURN
C
C        CODE FOR INCREMENTS EQUAL TO 1.
C
   20 XMAX = ABS(X(1))
      DO 30 I = 2,N
          XMAG = ABS(X(I))
          IF(XMAG.LE.XMAX) GO TO 30
          IDAMAX = I
          XMAX = XMAG
   30 CONTINUE
      RETURN
      END