text:
      subroutine bc3LHS (iBC,  BC,  iens,  EGmass )
c
c----------------------------------------------------------------------
c
c This routine satisfies the BC of LHS mass matrix for a single
c element.
c
c input:
c  iBC   (nshg) 	: boundary condition code
c  BC    (nshg,11)      : Dirichlet BC constraint parameters
c  ien   (npro,nshl)	: ien array for this element
c  EGmass(npro,nedof,nedof) : element consistent mass matrix before BC
c
c output:
c  EGmass(npro,nedof,nedof): LHS mass matrix after BC is satisfied
c
c
c Farzin Shakib, Winter 1987.
c Zdenek Johan,  Spring 1990. (Modified for general divariant gas)
c----------------------------------------------------------------------
c
        include "common.h"
c
	dimension iBC(nshg),      ien(npro,nshl),
     &		  BC(nshg,11),    EGmass(npro,nedof,nedof)
	integer iens(npro,nshl)
c
c prefer to show explicit absolute value needed for cubic modes and
c higher rather than inline abs on pointer as in past versions
c iens is the signed ien array ien is unsigned
c
	ien=abs(iens)
c
c.... loop over elements
c
        do iel = 1, npro
c
c.... loop over number of shape functions for this element
c
	do inod = 1, nshl
c
c.... set up parameters
c
        in  = ien(iel,inod)
        if (iBC(in) .eq. 0) goto 5000
        ioff = (inod - 1) * nflow
        i1 = ioff + 1
        i2 = ioff + 2
        i3 = ioff + 3
        i4 = ioff + 4
        i5 = ioff + 5
c
c.... pressure
c
	if ( btest(iBC(in),2) ) then

           do i = 1, nedof
              EGmass(iel,i,i1) = zero
              EGmass(iel,i1,i) = zero
           enddo
           EGmass(iel,i1,i1) = one
	endif
c
c.... density  ( not activated yet for LHS matrix )
c

c
c.... velocities
c
c
c.... x1-velocity
c
        if ( ibits(iBC(in),3,3) .eq. 1) then
           do i = 1, nedof
              EGmass(iel,i3,i) = EGmass(iel,i3,i)
     &                         - BC(in,4) * EGmass(iel,i2,i)
              EGmass(iel,i4,i) = EGmass(iel,i4,i)
     &                         - BC(in,5) * EGmass(iel,i2,i)
           enddo
           do i = 1, nedof
              EGmass(iel,i,i3) = EGmass(iel,i,i3)
     &                         - BC(in,4) * EGmass(iel,i,i2)
              EGmass(iel,i,i4) = EGmass(iel,i,i4)
     &                         - BC(in,5) * EGmass(iel,i,i2)
           enddo
           do i = 1, nedof
              EGmass(iel,i,i2) = zero
              EGmass(iel,i2,i) = zero
           enddo
           EGmass(iel,i2,i2) = one
        endif
c
c.... x2-velocity
c
        if ( ibits(iBC(in),3,3) .eq. 2 ) then
           do i = 1, nedof
              EGmass(iel,i2,i) = EGmass(iel,i2,i)
     &                         - BC(in,4) * EGmass(iel,i3,i)
              EGmass(iel,i4,i) = EGmass(iel,i4,i)
     &                         - BC(in,5) * EGmass(iel,i3,i)
           enddo
           do i = 1, nedof
              EGmass(iel,i,i2) = EGmass(iel,i,i2)
     &                         - BC(in,4) * EGmass(iel,i,i3)
              EGmass(iel,i,i4) = EGmass(iel,i,i4)
     &                         - BC(in,5) * EGmass(iel,i,i3)
           enddo

           do i = 1, nedof
              EGmass(iel,i,i3) = zero
              EGmass(iel,i3,i) = zero
           enddo
           EGmass(iel,i3,i3) = one
        endif
c
c.... x1-velocity and x2-velocity
c
        if ( ibits(iBC(in),3,3) .eq. 3 ) then
           do i = 1, nedof
              EGmass(iel,i4,i) = EGmass(iel,i4,i)
     &                     - BC(in,4) * EGmass(iel,i2,i)
     &                     - BC(in,6) * EGmass(iel,i3,i)
           enddo
           do i = 1, nedof
              EGmass(iel,i,i4) = EGmass(iel,i,i4)
     &                     - BC(in,4) * EGmass(iel,i,i2)
     &                     - BC(in,6) * EGmass(iel,i,i3)
           enddo

           do i = 1, nedof
              EGmass(iel,i,i2) = zero
              EGmass(iel,i2,i) = zero
              EGmass(iel,i,i3) = zero
              EGmass(iel,i3,i) = zero
           enddo
           EGmass(iel,i2,i2) = one
           EGmass(iel,i3,i3) = one
        endif
c
c.... x3-velocity
c
        if ( ibits(iBC(in),3,3) .eq. 4 ) then
           do i = 1, nedof
              EGmass(iel,i2,i) = EGmass(iel,i2,i)
     &                         - BC(in,4) * EGmass(iel,i4,i)
              EGmass(iel,i3,i) = EGmass(iel,i3,i)
     &                         - BC(in,5) * EGmass(iel,i4,i)
           enddo
           do i = 1, nedof
              EGmass(iel,i,i2) = EGmass(iel,i,i2)
     &                         - BC(in,4) * EGmass(iel,i,i4)
              EGmass(iel,i,i3) = EGmass(iel,i,i3)
     &                         - BC(in,5) * EGmass(iel,i,i4)
           enddo

           do i = 1, nedof
              EGmass(iel,i,i4) = zero
              EGmass(iel,i4,i) = zero
           enddo
           EGmass(iel,i4,i4) = one
        endif
c
c.... x1-velocity and x3-velocity
c
        if ( ibits(iBC(in),3,3) .eq. 5 ) then
           do i = 1, nedof
              EGmass(iel,i3,i) = EGmass(iel,i3,i)
     &                     - BC(in,4) * EGmass(iel,i2,i)
     &                     - BC(in,6) * EGmass(iel,i4,i)
           enddo
           do i = 1, nedof
              EGmass(iel,i,i3) = EGmass(iel,i,i3)
     &                     - BC(in,4) * EGmass(iel,i,i2)
     &                     - BC(in,6) * EGmass(iel,i,i4)
           enddo

           do i = 1, nedof
              EGmass(iel,i ,i2) = zero
              EGmass(iel,i2,i ) = zero
              EGmass(iel,i ,i4) = zero
              EGmass(iel,i4,i ) = zero
           enddo
           EGmass(iel,i2,i2) = one
           EGmass(iel,i4,i4) = one
        endif
c
c.... x2-velocity and x3-velocity
c
        if ( ibits(iBC(in),3,3) .eq. 6 ) then
           do i = 1, nedof
              EGmass(iel,i2,i) = EGmass(iel,i2,i)
     &                     - BC(in,4) * EGmass(iel,i3,i)
     &                     - BC(in,6) * EGmass(iel,i4,i)
           enddo
           do i = 1, nedof
              EGmass(iel,i,i2) = EGmass(iel,i,i2)
     &                     - BC(in,4) * EGmass(iel,i,i3)
     &                     - BC(in,6) * EGmass(iel,i,i4)
           enddo

           do i = 1, nedof
              EGmass(iel,i ,i3) = zero
              EGmass(iel,i3,i ) = zero
              EGmass(iel,i ,i4) = zero
              EGmass(iel,i4,i ) = zero
           enddo
           EGmass(iel,i3,i3) = one
           EGmass(iel,i4,i4) = one
        endif
c
c.... x1-velocity, x2-velocity, and x3-velocity
c
        if ( ibits(iBC(in),3,3) .eq. 7 ) then
           do i = 1, nedof
              EGmass(iel,i ,i2) = zero
              EGmass(iel,i2,i ) = zero
              EGmass(iel,i ,i3) = zero
              EGmass(iel,i3,i ) = zero
              EGmass(iel,i ,i4) = zero
              EGmass(iel,i4,i ) = zero
           enddo
           EGmass(iel,i2,i2) = one
           EGmass(iel,i3,i3) = one
           EGmass(iel,i4,i4) = one
        endif
c
c.... temperature
c
	if ( btest(iBC(in),1) ) then
           do i = 1, nedof
              EGmass(iel,i,i5) = zero
              EGmass(iel,i5,i) = zero
           enddo
           EGmass(iel,i5,i5) = one
	endif
c
c	Elaine
c.... scaled plane extraction boundary conditions
c
        if(intpres.eq.1) then
	  if ( btest(iBC(in),11) ) then
            do i = 1, nedof
              EGmass(iel,i ,i1) = zero
              EGmass(iel,i1,i ) = zero
              EGmass(iel,i ,i2) = zero
              EGmass(iel,i2,i ) = zero
              EGmass(iel,i ,i3) = zero
              EGmass(iel,i3,i ) = zero
              EGmass(iel,i ,i4) = zero
              EGmass(iel,i4,i ) = zero
              EGmass(iel,i,i5) = zero
              EGmass(iel,i5,i) = zero
            enddo
            EGmass(iel,i1,i1) = one
            EGmass(iel,i2,i2) = one
            EGmass(iel,i3,i3) = one
            EGmass(iel,i4,i4) = one
            EGmass(iel,i5,i5) = one
          endif
	else
	  if ( btest(iBC(in),11) ) then
            do i = 1, nedof
              EGmass(iel,i ,i2) = zero
              EGmass(iel,i2,i ) = zero
              EGmass(iel,i ,i3) = zero
              EGmass(iel,i3,i ) = zero
              EGmass(iel,i ,i4) = zero
              EGmass(iel,i4,i ) = zero
              EGmass(iel,i,i5) = zero
              EGmass(iel,i5,i) = zero
            enddo
            EGmass(iel,i2,i2) = one
            EGmass(iel,i3,i3) = one
            EGmass(iel,i4,i4) = one
            EGmass(iel,i5,i5) = one
          endif
	endif

 5000 continue

c
c.... end loop over shape functions (nodes)
c
      enddo
c
c.... end loop over elements
c
      enddo
c
c.... return
c
      return
      end
c
c
c
      subroutine bc3LHSSclr (iBC,  iens,  EGmasst )
c
c----------------------------------------------------------------------
c
c This routine satisfies the BC of LHS mass matrix for a single
c element.
c
c input:
c  iBC   (nshg) 	: boundary condition code
c  BC    (nshg,11)      : Dirichlet BC constraint parameters
c  ien   (npro,nshl)	: ien array for this element
c  EGmasst(npro,nshape,nshape) : element consistent mass matrix before BC
c
c output:
c  EGmasst(npro,nshape,nshape): LHS mass matrix after BC is satisfied
c
c
c Farzin Shakib, Winter 1987.
c Zdenek Johan,  Spring 1990. (Modified for general divariant gas)
c----------------------------------------------------------------------
c
        include "common.h"
c
	dimension iBC(nshg),      ien(npro,nshl),
     &		  EGmasst(npro,nshape,nshape)
c
	integer iens(npro,nshl)
c
c prefer to show explicit absolute value needed for cubic modes and
c higher rather than inline abs on pointer as in past versions
c iens is the signed ien array ien is unsigned
c
	ien=abs(iens)
c
        id = isclr+5
c.... loop over elements
c
        do iel = 1, npro
c
c.... loop over number of shape functions for this element
c
	do inod = 1, nshl
c
c.... set up parameters
c
        in  = ien(iel,inod)
        if (iBC(in) .eq. 0) goto 5000

c
c.... scalar
c
	if ( btest(iBC(in),id) ) then

           do i = 1, nshl
              EGmasst(iel,i,inod) = zero
              EGmasst(iel,inod,i) = zero
           enddo
           EGmasst(iel,inod,inod) = one
	endif

c

 5000 continue

c
c.... end loop over shape functions (nodes)
c
      enddo
c
c.... end loop over elements
c
      enddo
c
c.... return
c
      return
      end


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2021 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine beammpc(ipompc,nodempc,
     &  labmpc,nmpc,nmpc_,mpcfree,ikmpc,ilmpc,nk,nk_,nodeboun,ndirboun,
     &  ikboun,ilboun,nboun,nboun_,inode,node,co,typeboun)
!
!     generates an MPC for two nodes staying on a constant distance
!     from each other. The parameter inode indicates how many
!     times the present routine was called within the same *MPC
!     definition. For inode=1 "node" is node a, for inode=2 "node"
!     is node b.
!
      implicit none
!
      character*1 typeboun(*)
      character*20 labmpc(*)
!
      integer ipompc(*),nodempc(3,*),nmpc,nmpc_,mpcfree,nk,nk_,ikmpc(*),
     &  ilmpc(*),node,id,mpcfreeold,j,idof,l,nodeboun(*),nodea,nodeb,
     &  ndirboun(*),ikboun(*),ilboun(*),nboun,nboun_,inode,jmax,k,
     &  m
!
      real*8 co(3,*),dd,dmax
!
      save nodea
!
      if(inode.eq.1) then
         nodea=node
         return
      elseif(inode.eq.2) then
         nodeb=node
         dmax=0.d0
         do j=1,3
            dd=dabs(co(j,nodea)-co(j,nodeb))
            if(dd.gt.dmax) then
               dmax=dd
               jmax=j
            endif
         enddo
      endif
!
      nk=nk+1
      if(nk.gt.nk_) then
         write(*,*) '*ERROR in beammpc: increase nk_'
         call exit(201)
      endif
!
      j=jmax
      k=j+1
      if(k.gt.3) k=1
      l=k+1
      if(l.gt.3) l=1
!
      idof=8*(nodea-1)+j
      call nident(ikmpc,idof,nmpc,id)
      if(id.gt.0) then
         if(ikmpc(id).eq.idof) then
            write(*,*) '*WARNING in beammpc: DOF for node ',node
            write(*,*) '         in direction ',j,' has been used'
            write(*,*) '         on the dependent side of another MPC'
            write(*,*) '         PLANE constraint cannot be applied'
            return
         endif
      endif
      nmpc=nmpc+1
      if(nmpc.gt.nmpc_) then
         write(*,*) '*ERROR in beammpc: increase nmpc_'
         call exit(201)
      endif
!
      ipompc(nmpc)=mpcfree
      labmpc(nmpc)='BEAM                '
!
      do m=nmpc,id+2,-1
         ikmpc(m)=ikmpc(m-1)
         ilmpc(m)=ilmpc(m-1)
      enddo
      ikmpc(id+1)=idof
      ilmpc(id+1)=nmpc
!
      nodempc(1,mpcfree)=nodea
      nodempc(2,mpcfree)=j
      mpcfree=nodempc(3,mpcfree)
      nodempc(1,mpcfree)=nodea
      nodempc(2,mpcfree)=k
      mpcfree=nodempc(3,mpcfree)
      nodempc(1,mpcfree)=nodea
      nodempc(2,mpcfree)=l
      mpcfree=nodempc(3,mpcfree)
      nodempc(1,mpcfree)=nodeb
      nodempc(2,mpcfree)=j
      mpcfree=nodempc(3,mpcfree)
      nodempc(1,mpcfree)=nodeb
      nodempc(2,mpcfree)=k
      mpcfree=nodempc(3,mpcfree)
      nodempc(1,mpcfree)=nodeb
      nodempc(2,mpcfree)=l
      mpcfree=nodempc(3,mpcfree)
      nodempc(1,mpcfree)=nk
      nodempc(2,mpcfree)=j
      mpcfreeold=mpcfree
      mpcfree=nodempc(3,mpcfree)
      nodempc(3,mpcfreeold)=0
      idof=8*(nk-1)+j
      call nident(ikboun,idof,nboun,id)
      nboun=nboun+1
      if(nboun.gt.nboun_) then
         write(*,*) '*ERROR in beammpc: increase nboun_'
         call exit(201)
      endif
      nodeboun(nboun)=nk
      ndirboun(nboun)=j
      typeboun(nboun)='U'
      do m=nboun,id+2,-1
         ikboun(m)=ikboun(m-1)
         ilboun(m)=ilboun(m-1)
      enddo
      ikboun(id+1)=idof
      ilboun(id+1)=nboun
!
      return
      end

ext: f

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text: #redirect UC Davis Chancellor

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


PROGRAM target

  INTEGER, TARGET :: i, i2
  INTEGER, POINTER :: p_i, p_i2

  REAL, DIMENSION(:), ALLOCATABLE, TARGET :: rdata
  REAL, DIMENSION(:), POINTER :: pdata
  REAL, POINTER :: p_r

  PRINT *,'**** Simple variable ****'
  PRINT *,'Setting i to 10'
  i = 10
  p_i => i
  PRINT *, 'Pointer to i is now ', p_i
  PRINT *, 'Setting pointer to i to 20'
  p_i = 20
  PRINT *,'i is now ', i

  PRINT *, ''
  PRINT *, '**** Accidental equality ****'
  PRINT *, 'Setting i to 10'
  i = 10
  PRINT *, 'Setting i2 to 20'
  i2 = 20
  PRINT *, 'Pointing pointer to i'
  p_i => i
  PRINT *, 'Is pointer pointing to i :', ASSOCIATED(p_i, TARGET = i)
  PRINT *, 'Now using equality rather than points to operator'
  p_i = i2
  PRINT *, 'Is pointer still pointing to i :', ASSOCIATED(p_i, TARGET = i)
  PRINT *, 'p_i is ', p_i
  PRINT *, 'i is ', i
  PRINT *, 'i2 is ', i2
  PRINT *, 'p_i appears to correctly have the value 20 but by changing the &
      &value of i rather than pointing it to i2'
  PRINT *, 'Watch out for using = rather than => in pointer code'
  PRINT *, 'If you use = on a NULL pointer you will probably crash your code'
  PRINT *, 'Now setting i2 to 100'
  i2 = 100
  PRINT *, 'p_i is ', p_i
  PRINT *, 'i is ', i
  PRINT *, 'i2 is ', i2
  PRINT *, 'The value of p_i is not changed by changing i2'

  PRINT *,''
  PRINT *,'**** Arrays ****'
  PRINT *,'Allocating data to be -10:10'
  ALLOCATE(rdata(-10:10))
  pdata => rdata
  PRINT *, 'Bounds and sizes of data'
  PRINT *, LBOUND(rdata), UBOUND(rdata), SIZE(rdata)
  PRINT *, 'Bounds and sizes of pointer to data'
  PRINT *, LBOUND(pdata), UBOUND(pdata), SIZE(pdata)

  PRINT *, ''
  PRINT *, '**** Direct allocation ****'
  PRINT *,'Allocating p_i and setting to 100'
  ALLOCATE(p_i)
  p_i = 100

  PRINT *,'i is ', i
  PRINT *,'p_i is ', p_i
  PRINT *,'p_i is now just a variable and is no longer pointing to i'
  DEALLOCATE(p_i, rdata)

  PRINT *,''
  PRINT *,'**** Pointer to element of array ****'
  PRINT *,'Allocating array with 5 elements and setting to 0'
  ALLOCATE(rdata(1:5))
  rdata = 0
  PRINT *,'Pointing single pointer to element 3 and setting to 1'
  p_r => rdata(3)
  p_r = 1.0
  PRINT *,'Array contents is now'
  PRINT *, rdata
  DEALLOCATE(rdata)

  PRINT *,''
  PRINT *,'**** Pointer to array subsection ****'
  PRINT *,'Allocating array with 10 elements and setting to 0'
  ALLOCATE(rdata(1:10))
  rdata = 0
  PRINT *,'Setting array pointer to elements 3 to 5'
  pdata => rdata(3:5)
  PRINT '(A,I2,A,I2)','Size of pdata is ', LBOUND(pdata), ':', UBOUND(pdata)
  PRINT *,'Setting elements to 1'
  pdata = 1.0
  PRINT *,'Array contents is now'
  PRINT *, rdata
  DEALLOCATE(rdata)

  PRINT *, ''
  PRINT *, '**** No pointer to pointer ****'
  PRINT *, 'Setting i to 10 and i2 to 20'
  i  = 10
  i2 = 20
  PRINT *, 'Pointing p_i to i and p_i2 to p_i'
  p_i => i
  p_i2 => p_i

  PRINT *, 'Values of p_i and p_i2'
  PRINT *, 'p_i :', p_i
  PRINT *, 'p_i2 :', p_i2
  PRINT *, 'Pointing p_i to p2'
  p_i => i2
  PRINT *, 'Values of p_i and p_i2'
  PRINT *, 'p_i :', p_i
  PRINT *, 'p_i2 :', p_i2
  PRINT *, 'Is p_i2 pointing to i :', ASSOCIATED(p_i2, TARGET = i)
  PRINT *, 'Even though p_i2 was pointed to p_i it is NOT a pointer to p_i.'
  PRINT *, 'It is still a pointer to i'
  PRINT *, 'There is no concept of a pointer to a pointer in Fortran'

END PROGRAM target


text:
*----------------------------------------------------------------------*
      subroutine reorder_supvtx(possible,
     &     modify_contr,set_reord_list,reo_before,reo_info,
     &     contr,occ_vtx,rstr_vtx,idxop12,orb_info)
*----------------------------------------------------------------------*
*     check whether a single vertex is contracted to more than one
*     component of a super vertex
*     if this is the case, we may shift the involved occupations to
*     the super-vertex closest to the vertex they are commonly
*     contracted to, and we may add up the contraction
*
*     basically, this process is of relevance only if both contractions
*     refer to indices of the same C/A H/P/V/X (corresponds to a partial
*     symmetrization of the intermediate represented by the super
*     vertex); it will reduce the amount of storage necessary for the
*     intermediate
*----------------------------------------------------------------------*

      implicit none

      include 'opdim.h'
      include 'stdunit.h'
      include 'def_contraction.h'
      include 'def_reorder_info.h'
      include 'def_orbinf.h'
      include 'ifc_operators.h'

      integer, parameter ::
     &     ntest = 00

      logical, intent(out) ::
     &     possible
      logical, intent(in) ::
     &     modify_contr, set_reord_list, reo_before
      integer, intent(in) ::
     &     idxop12
      type(contraction), intent(inout) ::
     &     contr
      type(reorder_info), intent(inout) ::
     &     reo_info
      type(orbinf), intent(in), target ::
     &     orb_info
      integer, intent(inout) ::
     &     occ_vtx(ngastp,2,contr%nvtx)
      integer, intent(inout) ::
     &     rstr_vtx(2,orb_info%ngas,2,2,contr%nvtx)

      type(cntr_arc), pointer ::
     &     arc(:), arc_scr(:)
      type(cntr_vtx), pointer ::
     &     vertex(:)
      integer, pointer ::
     &     svertex(:), hpvxgas(:,:)
      integer ::
     &     occ_shr(ngastp,2), occ_shl(ngastp,2),
     &     cnt_shr(ngastp,2), cnt_shl(ngastp,2),
     &     rstr_shr(2,orb_info%ngas,2,2), rstr_shl(2,orb_info%ngas,2,2)
      integer ::
     &     narc, iarc, jarc, ivtx0, ivtx1, ivtx2, iprim, isuper,
     &     iarc_shift, iarc_sum, narc_new,
     &     maxreo, idx, idxsuper, ica, ica_vtx, hpvx,
     &     iblk, ivtx, idxnew, ireo, ngas, nspin
      logical ::
     &     renamed(contr%nvtx), modified
      logical, pointer ::
     &     reo_generated(:)
      integer, pointer ::
     &     arc_involved(:)

      integer, external ::
     &     imltlist

c dbg
c      print *,'reorder_supvtx: on input nreo ',reo_info%nreo
c dbg
      if (ntest.ge.100) then
        call write_title(lulog,wst_dbg_subr,'reorder_supvtx')
        write(lulog,*) 'modify_contr: ',modify_contr
        write(lulog,*) 'set_reord_list: ',set_reord_list
        write(lulog,*) 'idxop12: ',idxop12
        write(lulog,*) 'contr on input:'
        call prt_contr3(lulog,contr,occ_vtx)
      end if

      if (orb_info%nspin.gt.1)
     &     call quit(1,'reorder_supvtx','adapt for nspin>1')
      nspin = orb_info%nspin
      ngas = orb_info%ngas
      hpvxgas => orb_info%ihpvgas

      possible = .true.

      maxreo = 2*contr%nvtx
      if (set_reord_list.and..not.associated(reo_info%reo)) then
        reo_info%nreo = 0
        allocate(reo_info%reo(maxreo),reo_info%nca_vtx(contr%nvtx))
        reo_info%nvtx_contr = contr%nvtx
        call set_nca_vtx(reo_info%nca_vtx,occ_vtx,contr%nvtx)
      end if
      allocate(reo_generated(maxreo),arc_involved(maxreo))
      reo_generated = .false.
      arc_involved = -1

      narc = contr%narc
      if (modify_contr) then
        arc => contr%arc
      else
        allocate(arc_scr(contr%narc))
        arc_scr = contr%arc
        arc => arc_scr
      end if
      svertex => contr%svertex
      vertex =>  contr%vertex
      ! loop over arcs
      do iarc = 1, narc
        ! deleted arc?
        if (arc(iarc)%link(1).le.0) cycle
        ! loop over arcs and look for identical link(1) (if jarc<iarc)
        !                      or for identical link(2) (if jarc>iarc)
        do jarc = 1, narc
          ! deleted arc? or the arc itself?
          if (arc(iarc)%link(1).le.0.or.iarc.eq.jarc) cycle

c dbg
c          print *,'iarc, jarc: ',iarc,jarc
c dbg

          ! primitive vertex
          if (iarc.gt.jarc) iprim = 1
          if (iarc.lt.jarc) iprim = 2
          ! primitive vertex is the same?
c dbg
c          print *,'iprim = ',iprim
c          print *,'arc(iarc)%link(iprim),arc(jarc)%link(iprim):',
c     &         arc(iarc)%link(iprim),arc(jarc)%link(iprim)
c dbg
          if (arc(iarc)%link(iprim).ne.arc(jarc)%link(iprim)) cycle

          ! assumed super vertex
          if (iarc.gt.jarc) isuper = 2
          if (iarc.lt.jarc) isuper = 1

          ivtx0 = arc(iarc)%link(iprim)
          ivtx1 = arc(iarc)%link(isuper)
          ivtx2 = arc(jarc)%link(isuper)

c fix -- do not consider this case:
c          if (vertex(ivtx1)%idx_op.ne.idxop12) cycle

c dbg
c          print *,'iprim,isuper:',iprim,isuper
c          print *,'ivtx0,1,2: ',ivtx0,ivtx1,ivtx2
c dbg
          ! ensure ivtx0 < ivtx1 and ivtx0 < ivtx2
          !     or ivtx0 > ivtx1 and ivtx0 > ivtx2
          if (.not. ( (ivtx0.lt.ivtx1.and.ivtx0.lt.ivtx2) .or.
     &                (ivtx0.gt.ivtx1.and.ivtx0.gt.ivtx2) ) ) cycle
c dbg
c          print *,'survived a)'
c dbg

          ! same super vertex?
          if (svertex(ivtx1).ne.svertex(ivtx2)) cycle
c dbg
c          print *,'survived b)'
c dbg

          ! do contractions have non-vanishing overlap?
c          if (iocc_zero(iocc_overlap(arc(iarc)%occ_cnt,.false.,
c     &                               arc(jarc)%occ_cnt,.false.))) cycle

          occ_shr = 0
          occ_shl = 0
          cnt_shr = 0
          cnt_shl = 0
          do ica = 1, 2
            ! which C/A shall we look at in the associated vertices?
            ica_vtx = ica
            if (ivtx0.lt.ivtx1) ica_vtx = 3-ica
            do hpvx = 1, ngastp
              if (arc(iarc)%occ_cnt(hpvx,ica).gt.0 .and.
     &            arc(jarc)%occ_cnt(hpvx,ica).gt.0) then
                if (arc(iarc)%occ_cnt(hpvx,ica).eq.
     &                        occ_vtx(hpvx,ica_vtx,ivtx1)) then
                  ! iarc fully contracted to ivtx1
                  ! -> move jarc here
c dbg
c                  print *,'left case'
c dbg

                  cnt_shl(hpvx,ica) = arc(jarc)%occ_cnt(hpvx,ica)
                  occ_shl(hpvx,ica_vtx) = arc(jarc)%occ_cnt(hpvx,ica)

                else if (arc(jarc)%occ_cnt(hpvx,ica).eq.
     &                        occ_vtx(hpvx,ica_vtx,ivtx2)) then
                  ! jarc fully contracted to ivtx2
                  ! -> move iarc here
c dbg
c                  print *,'right case'
c dbg

                  cnt_shr(hpvx,ica) = arc(iarc)%occ_cnt(hpvx,ica)
                  occ_shr(hpvx,ica_vtx) = arc(iarc)%occ_cnt(hpvx,ica)

                else if (contr%nsupvtx.eq.2) then
                ! former treatment (outcommented by cmh) lead to problems
                ! in case where 2 njoined=2 intermediates were involved
                ! (leading to a njoined=2 result)
c dbg
c                  print *,'special case: just going ahead'
c dbg

cmh                  cnt_shr(hpvx,ica) = arc(iarc)%occ_cnt(hpvx,ica)
cmh                  occ_shr(hpvx,ica_vtx) = arc(iarc)%occ_cnt(hpvx,ica)

cmh                  occ_shr(hpvx,ica_vtx) = occ_vtx(hpvx,ica_vtx,ivtx1)
cmh     &                 - arc(iarc)%occ_cnt(hpvx,ica)

                else
c dbg
c                  print *,'skipping difficult reo'
c dbg
cmh                  possible = .false.
cmh               go ahead for now, maybe it will become possible after
cmh               some other contractions. If it remains not allowed,
cmh               this will be checked by allowed_contr
cmh                  return
c                  call quit(1,'reorder_supvtx','not yet')

                end if
              end if

            end do
          end do
c dbg
c          print *,'CNT SHL:'
c          call wrt_occ(lulog,cnt_shl)
c          print *,'CNT SHR:'
c          call wrt_occ(lulog,cnt_shr)
c          print *,'OCC SHL:'
c          call wrt_occ(lulog,occ_shl)
c          print *,'OCC SHR:'
c          call wrt_occ(lulog,occ_shr)
c dbg

          ! update contractions
          arc(iarc)%occ_cnt =
     &         arc(iarc)%occ_cnt + cnt_shl - cnt_shr
          arc(jarc)%occ_cnt =
     &         arc(jarc)%occ_cnt - cnt_shl + cnt_shr

          ! update super-vertex occupation
          if (modify_contr) then
            call fit_restr(rstr_shr,occ_shr,
     &           rstr_vtx(1,1,1,1,ivtx1),hpvxgas,ngas)
            call fit_restr(rstr_shl,occ_shl,
     &           rstr_vtx(1,1,1,1,ivtx2),hpvxgas,ngas)
c dbg
c            print *,'vtx # 1:'
c            call wrt_occ_rstr(6,ivtx1,occ_vtx(1,1,ivtx1),
c     &                        rstr_vtx(1,1,1,1,ivtx1),ngas,nspin)
c            print *,'vtx # 2:'
c            call wrt_occ_rstr(6,ivtx2,occ_vtx(1,1,ivtx2),
c     &                        rstr_vtx(1,1,1,1,ivtx2),ngas,nspin)
c            print *,'occ_shr'
c            call wrt_occ_rstr(6,0,occ_shr,
c     &                        rstr_shr,ngas,nspin)
c            print *,'occ_shl'
c            call wrt_occ_rstr(6,0,occ_shl,
c     &                        rstr_shl,ngas,nspin)
c dbg
c            rstr_vtx(1:2,1:ngas,1:2,1:2,1:nspin,ivtx1) =
c     &           rstr_vtx(1:2,1:ngas,1:2,1:2,1:nspin,ivtx1)
            rstr_vtx(1:2,1:ngas,1:2,1:2,ivtx1) =
     &           rstr_vtx(1:2,1:ngas,1:2,1:2,ivtx1)
     &           + rstr_shl - rstr_shr
c            rstr_vtx(1:2,1:ngas,1:2,1:2,1:nspin,ivtx2) =
c     &           rstr_vtx(1:2,1:ngas,1:2,1:2,1:nspin,ivtx2)
            rstr_vtx(1:2,1:ngas,1:2,1:2,ivtx2) =
     &           rstr_vtx(1:2,1:ngas,1:2,1:2,ivtx2)
     &           - rstr_shl + rstr_shr
            occ_vtx(1:ngastp,1:2,ivtx1) =
     &           occ_vtx(1:ngastp,1:2,ivtx1) + occ_shl - occ_shr
            occ_vtx(1:ngastp,1:2,ivtx2) =
     &           occ_vtx(1:ngastp,1:2,ivtx2) - occ_shl + occ_shr
c dbg
c            print *,'new vtx # 1:'
c            call wrt_occ_rstr(6,ivtx1,occ_vtx(1,1,ivtx1),
c     &                        rstr_vtx(1,1,1,1,ivtx1),ngas,nspin)
c            print *,'new vtx # 2:'
c            call wrt_occ_rstr(6,ivtx2,occ_vtx(1,1,ivtx2),
c     &                        rstr_vtx(1,1,1,1,ivtx2),ngas,nspin)
c dbg
          end if

          if (set_reord_list) then
            if (iocc_nonzero(occ_shr)) then
              reo_info%nreo = reo_info%nreo+1
              reo_generated(reo_info%nreo) = .true.
              arc_involved(reo_info%nreo) = iarc
              idx = reo_info%nreo
              if (idx.gt.maxreo) then
                write(lulog,*) 'idx,maxreo: ',idx,maxreo
                call quit(1,'reorder_supvtx','unexpected event')
              end if
              idxsuper = svertex(ivtx1)
              reo_info%reo(idx)%idxsuper = idxsuper
              reo_info%reo(idx)%idxop_ori =
     &             vertex(ivtx1)%idx_op
              reo_info%reo(idx)%iblkop_ori =
     &             vertex(ivtx1)%iblk_op
              reo_info%reo(idx)%dagger_ori =
     &             vertex(ivtx1)%dagger
              ! flag whether this is the result vertex of prev. binary contr.
              reo_info%reo(idx)%is_bc_result =
     &             idxop12.eq.vertex(ivtx1)%idx_op
              ! flag whether reo should occur before contraction
              reo_info%reo(idx)%reo_before   = reo_before
              ! which components of super-vertex
              reo_info%reo(idx)%to =
     &             imltlist(idxsuper,svertex,ivtx2,1)
              reo_info%reo(idx)%from =
     &             imltlist(idxsuper,svertex,ivtx1,1)
              reo_info%reo(idx)%to_vtx   = ivtx2
              reo_info%reo(idx)%from_vtx = ivtx1
              ! shifted occupation
              reo_info%reo(idx)%occ_shift = occ_shr
              reo_info%reo(idx)%shift_i0 = .false.
            end if
            if (iocc_nonzero(occ_shl)) then
              reo_info%nreo = reo_info%nreo+1
              reo_generated(reo_info%nreo) = .true.
              arc_involved(reo_info%nreo) = jarc
              idx = reo_info%nreo
              if (idx.gt.maxreo) then
                write(lulog,*) 'idx,maxreo: ',idx,maxreo
                call quit(1,'reorder_supvtx','unexpected event')
              end if
              idxsuper = svertex(ivtx1)
              reo_info%reo(idx)%idxsuper = idxsuper
              reo_info%reo(idx)%idxop_ori =
     &             vertex(ivtx1)%idx_op
              reo_info%reo(idx)%iblkop_ori =
     &             vertex(ivtx1)%iblk_op
              reo_info%reo(idx)%dagger_ori =
     &             vertex(ivtx1)%dagger
              ! flag whether this is the result vertex of prev. binary contr.
              reo_info%reo(idx)%is_bc_result =
     &             idxop12.eq.vertex(ivtx1)%idx_op
              ! flag whether reo should occur before contraction
              reo_info%reo(idx)%reo_before   = reo_before
              ! which components of super-vertex
              reo_info%reo(idx)%to =
     &             imltlist(idxsuper,svertex,ivtx1,1)
              reo_info%reo(idx)%from =
     &             imltlist(idxsuper,svertex,ivtx2,1)
              reo_info%reo(idx)%to_vtx   = ivtx1
              reo_info%reo(idx)%from_vtx = ivtx2
              ! shifted occupation
              reo_info%reo(idx)%occ_shift = occ_shl
              reo_info%reo(idx)%shift_i0 = .false.
            end if

          end if

        end do

      end do

      ! remove deleted or zero arcs:
      if (modify_contr) then
        modified = .false.
        narc_new = 0
        do iarc = 1, narc
          if (arc(iarc)%link(1).le.0 .or.
     &        iocc_zero(arc(iarc)%occ_cnt) ) cycle
          narc_new = narc_new+1
          if (narc_new.lt.iarc) arc(narc_new) = arc(iarc)
        end do
        contr%narc = narc_new
        modified = narc.ne.narc_new

        ! a quickie: new intermediate
        renamed = .false.
        idxnew = idxop12-1
        do ireo = 1, reo_info%nreo
          if (.not.reo_generated(ireo)) cycle

          modified = .true.

          reo_info%reo(ireo)%idxop_new = idxop12 ! default
          reo_info%reo(ireo)%dagger_new = .false. ! default
          idxsuper = reo_info%reo(idx)%idxsuper
          iblk = 0
          do ivtx = 1, contr%nvtx
            if (renamed(ivtx)) cycle ! do not rename twice
            if (svertex(ivtx).ne.idxsuper) cycle
            if (reo_before.or.vertex(ivtx)%idx_op.ne.idxop12) then
              renamed(ivtx) = .true.
              vertex(ivtx)%idx_op = idxnew
              iblk = iblk+1
              vertex(ivtx)%iblk_op = iblk
              vertex(ivtx)%dagger = .false.
            end if
          end do
          if (iblk.ne.0) reo_info%reo(ireo)%idxop_new = idxnew
          if (iblk.ne.0) idxnew = idxnew-1

          if (contr%index_info)
     &         call reorder_string_info(contr,
     &              reo_info%reo(ireo)%from_vtx,
     &              reo_info%reo(ireo)%to_vtx,
     &              arc_involved(ireo),
     &              reo_info%reo(ireo)%occ_shift,
     &              .false.)

        end do

        ! set 0-contraction, if necessary
        ! as they were removed above
        call check_disconnected(contr)

        if (contr%index_info.and.modified)
     &       call update_string_info(contr)

      end if

      if (.not.modify_contr) deallocate(arc_scr)

      deallocate(reo_generated,arc_involved)

      if (ntest.ge.100) then
        write(lulog,*) 'contr at the end of reorder_supvtx: nreo = ',
     &       reo_info%nreo
        if (.not.modify_contr) write(lulog,*) 'should not have changed!'
        call prt_contr3(lulog,contr,occ_vtx)
      end if

      end

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!  contact: delft3d.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands
!!
!!  All indications and logos of, and references to registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

      subroutine ddepth ( pmsa   , fl     , ipoint , increm , noseg  ,
     &                    noflux , iexpnt , iknmrk , noq1   , noq2   ,
     &                    noq3   , noq4   )
!>\file
!>       Dynamic calculation of the depth as volume / surf

!
!     Description of the module :
!
!        General water quality module for DELWAQ:
!        DEPTH CALCULATION FROM HORIZONTAL SURFACE AREA OF A SEGMENT
!
! Name    T   L I/O   Description                                    Units
! ----    --- -  -    -------------------                             ----
! DEPTH   R*4 1 O depth of the water column                            [m]
! SURF    R*4 1 I surface area of segment                             [m2]
! VOLUME  R*4 1 I volume of segment                                   [m3]

!     Logical Units : -

!     Modules called : -

!     Name     Type   Library
!     ------   -----  ------------

      IMPLICIT REAL (A-H,J-Z)

      REAL     PMSA  ( * ) , FL    (*)
      INTEGER  IPOINT( * ) , INCREM(*) , NOSEG , NOFLUX,
     +         IEXPNT(4,*) , IKNMRK(*) , NOQ1, NOQ2, NOQ3, NOQ4
      character(55) message

      message = 'SURF in DDEPTH zero at segment:'
      IP1  = IPOINT( 1)
      IP2  = IPOINT( 2)
      IP3  = IPOINT( 3)
!
      IFLUX = 0
      DO 9000 ISEG = 1 , NOSEG
      CALL DHKMRK(3,IKNMRK(ISEG),IKMRK3)
      IF (IKMRK3.EQ.1 .OR. IKMRK3.EQ.3) THEN
!
      VOLUME = PMSA(IP1 )
      SURF   = PMSA(IP2 )

      IF (SURF    .LT. 1E-30) THEN
         write ( message(32:55) , '(i9,1x,e14.6)' ) iseg, surf
         CALL ERRSYS ( message, 1 )
      ENDIF

!***********************************************************************
!**** Calculate DEPTH
!***********************************************************************
!
      DEPTH = VOLUME / SURF
!
      PMSA (IP3 ) = DEPTH
!
      ENDIF
!
      IFLUX = IFLUX + NOFLUX
      IP1   = IP1   + INCREM (  1 )
      IP2   = IP2   + INCREM (  2 )
      IP3   = IP3   + INCREM (  3 )
!
 9000 CONTINUE
!
      RETURN
!
      END

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
Describe Users/LindsayHartman here.

20080227 16:30:12 nbsp Welcome to the Wiki. After reading your http://daviswiki.org/Chiles_Mansion?actiondiff&version226&version125 comment on the Chiles Mansion page Id like to invite you to add more information about the Native Peoples who were here before. Remember that the Wiki is only as good as you make it. Oh, and generally comments left in the comment section should have attribution for who left them to make it easier to track the history of the conversation; to do this you can use the third button from the right on the edit page (the one that looks like a name tag). Users/JasonAller


text:
c-----------------------------------------------------------------------
      subroutine proxy_setup(a,b,c,d,z,w,g)

      include 'SIZE'
      include 'TOTAL'

      real a(lx1*lx1),b(lx1),c(lx1*lx1),d(lx1*lx1),z(lx1)
     $               , w(lx1*2),g(nx1*ny1*nz1,6,nelt)

      call semhat(a,b,c,d,z,w,nx1-1)

      n = nx1*nx1
      call copy(dxm1,d,n)
      call transpose(dxtm1,nx1,dxm1,nx1)

      call copy(zgm1,z,nx1)   ! GLL points
      call copy(wxm1,b,nx1)   ! GLL weights

      call setup_g(g)

c     m = nx1*ny1*nz1*nelt
c     call outmat(g,6,m,'gxyz 1',m)

      return
      end
c-------------------------------------------------------------------------
      subroutine setup_g(g)

      include 'SIZE'
      include 'TOTAL'
      real g(nx1,ny1,nz1,6,nelt)
      integer e

      n = nx1*ny1*nz1*nelt

      g = 0

      do e=1,nelt
      do k=1,nz1
      do j=1,ny1
      do i=1,nx1
c        call rzero(g(1,i,j,k,e),6)
         g(i,j,k,1,e) = wxm1(i)*wxm1(j)*wxm1(k)
         g(i,j,k,4,e) = wxm1(i)*wxm1(j)*wxm1(k)
         g(i,j,k,6,e) = wxm1(i)*wxm1(j)*wxm1(k)
      enddo
      enddo
      enddo
      enddo

      return
      end
c-------------------------------------------------------------------------
      subroutine transpose(a,lda,b,ldb)
      real a(lda,1),b(ldb,1)
c
      do j=1,ldb
         do i=1,lda
            a(i,j) = b(j,i)
         enddo
      enddo
      return
      end
c-----------------------------------------------------------------------
      subroutine outmat(a,m,n,name6,ie)
      real a(m,n)
      character*6 name6
c
      n10 = min(n,10)
      write(6,*)
      write(6,*) ie,' matrix: ',name6,m,n
      do i=1,m
         write(6,6) ie,name6,(a(i,j),j=1,n10)
      enddo
    6 format(i3,1x,a6,1p10e12.4)
      write(6,*)
      return
      end
c-----------------------------------------------------------------------
      subroutine outmat1(a,m,n,name6,ie)
      real a(m,n)
      character*6 name6
c
      n10 = min(n,10)
      write(ie,*)
      write(ie,*) ie,' matrix: ',name6,m,n
      do i=1,m
         write(ie,6) ie,name6,(a(i,j),j=1,n10)
      enddo
    6 format(i3,1x,a6,1p10e12.4)
      write(ie,*)
      return
      end
c-----------------------------------------------------------------------
      function randx(seed)

      arg   = 1.e9*seed
      arg   = 1.e9*cos(arg)
      randx = sin(arg)
      seed  = randx
      seed  = randx

      return
      end
c-----------------------------------------------------------------------


text:
	SUBROUTINE GDPTDTM  ( gdatim, gvcord, gfunc, firstm, lastim,
     +			     time, iret )
C************************************************************************
C* GDPTDTM								*
C*									*
C* This subroutine gets the date/time requested by the user.		*
C*									*
C* GDPDTV  ( GDATIM, GVCORD, GFUNC, FIRSTM, LASTIM, TIME, IRET )	*
C*									*
C* Input parameters:							*
C*	GDATIM		CHAR*		User input for date/time	*
C*	GVCORD		CHAR*		User input for vert coord	*
C*	GFUNC		CHAR*		User input for grid function	*
C*	FIRSTM		CHAR*		First time in file		*
C*	LASTIM		CHAR*		Last time in file		*
C*									*
C* Output parameters:							*
C*	TIME  (2)	CHAR*		Date/time			*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					 -5 = invalid time		*
C**									*
C* Log:									*
C* M. Li/SAIC		08/07	Modified from GDPDTV			*
C************************************************************************
	CHARACTER*(*)	gdatim, gvcord, gfunc, firstm, lastim, time (2)
C*
	CHARACTER	ctime*48
C------------------------------------------------------------------------
	iret  = 0
C
C*	Read input grid function to see if it contains date/time.
C
	CALL ST_LSTR  ( gfunc, iend, ier )
	idt = INDEX ( gfunc, '^' )
	IF  ( idt .eq. 0 )  THEN
	    ctime = gdatim
	  ELSE
	    idt   = idt + 1
	    ict   = 0
	    ctime = ' '
	    DO WHILE  ( ( idt .le. iend ) .and.
     +			( gfunc (idt:idt) .ne. ')' ) .and.
     +			( gfunc (idt:idt) .ne. ']' ) .and.
     +			( gfunc (idt:idt) .ne. ',' ) .and.
     +			( gfunc (idt:idt) .ne. ';' ) .and.
     +			( gfunc (idt:idt) .ne. '^' ) .and.
     +			( gfunc (idt:idt) .ne. '%' ) .and.
     +			( gfunc (idt:idt) .ne. '+' ) .and.
     +			( gfunc (idt:idt) .ne. '@' ) )
		ict = ict + 1
		ctime (ict:ict) = gfunc (idt:idt)
		idt = idt + 1
	    END DO
	END IF
C
C*	Change time into two parts.
C
	CALL GR_GTIM  ( ctime, firstm, lastim, time (1), time (2), ier )
	IF  ( ier .ne. 0 )  THEN
	    iret = -5
	    RETURN
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
great_rsd(v1,kk1).
great_rsd(cc1,aa1).
great_rsd(a1,dd1).
great_rsd(f1,v1).
great_rsd(cc1,m1).
great_rsd(cc1,k1).
great_rsd(ff1,m1).
great_rsd(f1,dd1).
great_rsd(e1,l1).
great_rsd(u1,l1).
great_rsd(cc1,z1).
great_rsd(aa1,h1).
great_rsd(p1,z1).
great_rsd(cc1,l1).
great_rsd(v1,x1).
great_rsd(e1,p1).
great_rsd(q1,dd1).
great_rsd(e1,g1).
great_rsd(ff1,hh1).
great_rsd(e1,c1).
great_rsd(p1,w1).
great_rsd(f1,w1).
great_rsd(e1,hh1).
great_rsd(u1,g1).
great_rsd(ee1,jj1).
great_rsd(aa1,jj1).
great_rsd(t1,g1).
great_rsd(m1,h1).
great_rsd(z1,q1).
great_rsd(ii1,u1).
great_rsd(g1,e1).
great_rsd(ff1,aa1).

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


text:
module dim
      integer, parameter:: ndim=3
end module dim
module numz
! module defines the basic real types
      integer, parameter:: b8 = selected_real_kind(14)
      integer, parameter:: b16 = selected_real_kind(14)
      integer, parameter :: cin=5
      integer, parameter :: cout=6
      integer, parameter :: cerr=0
      character(len=11):: form1="(1x,e22.15)"
end module

module face
interface
    subroutine get_acc_jerk_pot_coll(mass, pos, &
                                     vel,  acc, &
                                     jerk,  n, epot, &
                                     coll_time)
    use numz
    implicit none
    real(b8) mass(0:), pos(0:,0:)
    real(b8) vel(0:,0:),  acc(0:,0:)
    real(b8) jerk(0:,0:)
    integer  n
    real(b8) epot
    real(b8) coll_time
    end subroutine

    subroutine evolve_step(mass,  pos,  vel, &
                  acc,  jerk, n, t, &
                  dt,  epot, coll_time)
	use numz
	use dim
	implicit none
    real(b8) mass(0:), pos(0:,0:)
    real(b8) vel(0:,0:),  acc(0:,0:)
    real(b8) jerk(0:,0:)
    integer  n
    real(b8) t,dt
    real(b8) epot
    real(b8) coll_time
    end subroutine

    subroutine predict_step(pos,vel, acc, jerk,n,dt)
	use numz
	use dim
	implicit none
    real(b8) pos(0:,0:)
    real(b8) vel(0:,0:)
    real(b8) acc(0:,0:)
    real(b8) jerk(0:,0:)
    integer  n
    real(b8) dt
    end subroutine

    subroutine correct_step(pos,vel, acc, jerk, &
                  old_pos, old_vel, old_acc, old_jerk, &
                  n,dt)
	use numz
	use dim
	implicit none
    real(b8) pos(0:,0:),vel(0:,0:),acc(0:,0:),jerk(0:,0:)
    real(b8) old_pos(0:,0:),old_vel(0:,0:),old_acc(0:,0:),old_jerk(0:,0:)
    integer  n
    real(b8) dt
    end subroutine

    subroutine put_snapshot(mass, pos, vel, n, t)
 	use numz
	use dim
	implicit none
    real(b8) mass(0:), pos(0:,0:)
    real(b8) vel(0:,0:)
    integer  n
    real(b8) t
    end subroutine

    subroutine write_diagnostics(mass, pos, &
                       vel, acc, &
                       jerk, n, t, epot, &
                       nsteps, einit,init_flag, &
                       x_flag)
	use numz
	use dim
	implicit none
    real(b8) mass(0:), pos(0:,0:)
    real(b8) vel(0:,0:),  acc(0:,0:)
    real(b8) jerk(0:,0:)
    integer  n
    real(b8) t
    real(b8) epot
    integer nsteps
    real(b8) einit
    logical init_flag,x_flag
    end subroutine

    subroutine read_options(flag, dt_param, dt_dia, dt_out, dt_tot, i_flag,x_flag)
	use numz
	use dim
	implicit none
    logical flag
    real(b8) dt_param,dt_dia,dt_out,dt_tot
    logical i_flag,x_flag
    end subroutine

    subroutine get_snapshot(mass,  pos,  vel, n)
	use numz
	use dim
	implicit none
    real(b8) mass(0:), pos(0:,0:)
    real(b8) vel(0:,0:)
    integer n
    end subroutine

    subroutine evolve(mass, pos, vel, &
                  n,  t, dt_param, dt_dia, dt_out, &
                  dt_tot,  init_out,  x_flag)
	use numz
	use dim
	implicit none
    real(b8) mass(0:), pos(0:,0:)
    real(b8) vel(0:,0:)
    integer  n
    real(b8) t, dt_param, dt_dia, dt_out, dt_tot
    logical init_out,x_flag
    end subroutine

end interface
end module face

!/*-----------------------------------------------------------------------------
! *  get_acc_jerk_pot_coll  --  calculates accelerations and jerks, and as side
! *                             effects also calculates potential energy and
! *                             the time scale coll_time for significant changes
! *                             in local configurations to occur.
! *                                                  __                     __
! *                                                 |          -->  -->       |
! *               M                           M     |           r  . v        |
! *   -->          j    -->       -->          j    | -->        ji   ji -->  |
! *    a   ==  --------  r    ;    j   ==  -------- |  v   - 3 ---------  r   |
! *     ji     |-->  |3   ji        ji     |-->  |3 |   ji      |-->  |2   ji |
! *            | r   |                     | r   |  |           | r   |       |
! *            |  ji |                     |  ji |  |__         |  ji |     __|
! *
! *  note: it would be cleaner to calculate potential energy and collision time
! *        in a separate function.  However, the current function is by far the
! *        most time consuming part of the whole program, with a double loop
! *        over all particles that is executed every time step.  Splitting off
! *        some of the work to another function would significantly increase
! *        the total computer time (by an amount close to a factor two).
! *
! *  We determine the values of all four quantities of interest by walking
! *  through the system in a double {i,j} loop.  The first three, acceleration,
! *  jerk, and potential energy, are calculated by adding successive terms;
! *  the last, the estimate for the collision time, is found by determining the
! *  minimum value over all particle pairs and over the two choices of collision
! *  time, position/velocity and sqrt(position/acceleration), where position and
! *  velocity indicate their relative values between the two particles, while
! *  acceleration indicates their pairwise acceleration.  At the start, the
! *  first three quantities are set to zero, to prepare for accumulation, while
! *  the last one is set to a very large number, to prepare for minimization.
! *       The integration loops only over half of the pairs, with j > i, since
! *  the contributions to the acceleration and jerk of particle j on particle i
! *  is the same as those of particle i on particle j, apart from a minus sign
! *  and a different mass factor.
! *-----------------------------------------------------------------------------
! */

subroutine get_acc_jerk_pot_coll(mass, pos, &
                                 vel,  acc, &
                                 jerk, n, epot, &
                                 coll_time)
    use numz
    use dim
    implicit none
    real(b8) mass(0:), pos(0:,0:)
    real(b8) vel(0:,0:),  acc(0:,0:)
    real(b8) jerk(0:,0:)
    integer  n
    real(b8) epot
    real(b8) coll_time
    real(b8) rji(0:ndim-1)
    real(b8) vji(0:ndim-1)
    real r2,v2,rv_r2,r,r3,da(0:2),dj(0:2),da2
    real(b8) coll_time_q ! collision time to 4th power
    real(b8) coll_est_q ! collision time scale estimate to 4th power (quartic)
    real(b8), parameter :: vln =1.0e300_b8
    real(b16) mij
    integer i,k,j

    acc=0.0_b8
    jerk=0.0_b8
    epot=0.0_b8
    coll_time_q = vln      ! collision time to 4th power
    do i=0,n-1
      do j=i+1, n-1                              !! rji() is the vector from
!            real rji(NDIM)                        !! particle i to particle j
!            real vji(NDIM)                        !! vji() = d rji() / d t
            do k=0, ndim-1
                rji(k) = pos(j,k) - pos(i,k)
                vji(k) = vel(j,k) - vel(i,k)
            enddo
            r2 = 0.0_b8                           !! | rji |^2
            v2 = 0.0_b8                           !! | vji |^2
            rv_r2 = 0.0_b8                        !! ( rij . vij ) / | rji |^2
            do k=0, ndim-1
                r2=r2+ rji(k) * rji(k)
                v2=v2+ vji(k) * vji(k)
                rv_r2= rv_r2+ rji(k) * vji(k)
            enddo
            rv_r2 = rv_r2 / r2
            r = sqrt(r2)                     !! | rji |
            r3 = r * r2                      !! | rji |^3

!! add the {i,j} contribution to the total potential energy for the system:

            epot = epot - mass(i) * mass(j) / r

!! add the {j (i)} contribution to the {i (j)} values of acceleration and jerk:

!           real da(3)                            !! main terms in pairwise
!           real dj(3)                            !! acceleration and jerk
            do k=0, ndim-1
                da(k) = rji(k) / r3                           !! see equations
                dj(k) = (vji(k) - 3 * rv_r2 * rji(k)) / r3    !! in the header
            enddo
            do k=0, ndim-1
                acc(i,k) = acc(i,k) + mass(j) * da(k)                 !! using symmetry
                acc(j,k) = acc(j,k) - mass(i) * da(k)                 !! find pairwise
                jerk(i,k) = jerk(i,k) + mass(j) * dj(k)               !! acceleration
                jerk(j,k) = jerk(j,k) - mass(i) * dj(k)               !! and jerk
            enddo

!! first collision time estimate, based on unaccelerated linear motion:

            coll_est_q = (r2*r2) / (v2*v2)
            if (coll_time_q .gt. coll_est_q) coll_time_q = coll_est_q

!! second collision time estimate, based on free fall:

            da2 = 0.0_b8                                  !! da2 becomes the
            do k=0, ndim-1                              !! square of the
                da2= da2+ da(k) * da(k)
            enddo                                         !! pair-wise accel-
            mij = mass(i) + mass(j)                       !! eration between
            da2=da2* mij * mij                            !! particles i and j

            coll_est_q = r2/da2
            if (coll_time_q .gt. coll_est_q) coll_time_q = coll_est_q
        enddo
    enddo                                              !! from q for quartic back
    coll_time = sqrt(sqrt(coll_time_q))                !! to linear collision time
end subroutine get_acc_jerk_pot_coll

subroutine evolve_step(mass,  pos,  vel, &
                  acc,  jerk, n, t, &
                  dt,  epot, coll_time)
	use numz
	use dim
	use face, only:predict_step,get_acc_jerk_pot_coll,correct_step
	implicit none
    real(b8) mass(0:), pos(0:,0:)
    real(b8) vel(0:,0:),  acc(0:,0:)
    real(b8) jerk(0:,0:)
    integer  n
    real(b8) t,dt
    real(b8) epot
    real(b8) coll_time

	real(b8),allocatable:: old_pos(:,:),old_vel(:,:),old_acc(:,:),old_jerk(:,:)
	allocate( old_pos(0:n-1,0:ndim-1) )
	allocate( old_vel(0:n-1,0:ndim-1) )
	allocate( old_acc(0:n-1,0:ndim-1) )
	allocate( old_jerk(0:n-1,0:ndim-1) )
	old_pos=pos
	old_vel=vel
	old_acc=acc
	old_jerk=jerk

	call predict_step(pos, vel, acc, jerk, n, dt)
	call get_acc_jerk_pot_coll(mass, pos, vel, acc, jerk, n, epot, coll_time)
	call correct_step(pos, vel, acc, jerk, old_pos, old_vel, old_acc, old_jerk, n, dt)
	t = t +dt

	deallocate( old_pos     )
	deallocate( old_vel     )
	deallocate( old_acc     )
	deallocate( old_jerk     )
end subroutine  evolve_step

!/*-----------------------------------------------------------------------------
! *  predict_step  --  takes the first approximation of one Hermite integration
! *                    step, advancing the positions and velocities through a
! *                    Taylor series development up to the order of the jerks.
! *-----------------------------------------------------------------------------
! */

subroutine predict_step(pos,vel, acc, jerk,n,dt)
	use numz
	use dim
	implicit none
    real(b8) pos(0:,0:)
    real(b8) vel(0:,0:)
    real(b8) acc(0:,0:)
    real(b8) jerk(0:,0:)
    integer  n
    real(b8) dt

    real(b8) dt2,dt3
    integer i,k
    dt2=dt*dt
    dt3=dt2*dt
    dt2=dt2/2.0_b8
    dt3=dt3/6.0_b8
    do i=0,n-1
      do k=0,ndim-1

            pos(i,k) = pos(i,k) + vel(i,k)*dt + acc(i,k)*dt2 &
                                                 + jerk(i,k)*dt3
            vel(i,k) = vel(i,k) + acc(i,k)*dt + jerk(i,k)*dt2
      enddo
    enddo
end subroutine predict_step
!/*-----------------------------------------------------------------------------
! *  correct_step  --  takes one iteration to improve the new values of position
! *                    and velocities, effectively by using a higher-order
! *                    Taylor series constructed from the terms up to jerk at
! *                    the beginning and the end of the time step.
! *-----------------------------------------------------------------------------
! */

subroutine correct_step(pos,vel, acc, jerk, &
                  old_pos, old_vel, old_acc, old_jerk, &
                  n,dt)
	use numz
	use dim
	implicit none
    real(b8) pos(0:,0:),vel(0:,0:),acc(0:,0:),jerk(0:,0:)
    real(b8) old_pos(0:,0:),old_vel(0:,0:),old_acc(0:,0:),old_jerk(0:,0:)
    integer  n
    real(b8) dt

    integer i,k

    do i=0,n-1
      do k=0,ndim-1
            vel(i,k) = old_vel(i,k) + (old_acc(i,k) + acc(i,k))*dt/2 &
                                      + (old_jerk(i,k) - jerk(i,k))*dt*dt/12
            pos(i,k) = old_pos(i,k) + (old_vel(i,k) + vel(i,k))*dt/2 &
                                      + (old_acc(i,k) - acc(i,k))*dt*dt/12
      enddo
    enddo
end subroutine correct_step
!/*-----------------------------------------------------------------------------
! *  evolve  --  integrates an N-body system, for a total duration dt_tot.
! *              Snapshots are sent to the standard output stream once every
! *              time interval dt_out.  Diagnostics are sent to the standard
! *              error stream once every time interval dt_dia.
! *
! *  note: the integration time step, shared by all particles at any given time,
! *        is variable.  Before each integration step we use coll_time (short
! *        for collision time, an estimate of the time scale for any significant
! *        change in configuration to happen), multiplying it by dt_param (the
! *        accuracy parameter governing the size of dt in units of coll_time),
! *        to obtain the new time step size.
! *
! *  Before moving any particles, we start with an initial diagnostics output
! *  and snapshot output if desired.  In order to write the diagnostics, we
! *  first have to calculate the potential energy, with get_acc_jerk_pot_coll().
! *  That function also calculates accelerations, jerks, and an estimate for the
! *  collision time scale, all of which are needed before we can enter the main
! *  integration loop below.
! *       In the main loop, we take as many integration time steps as needed to
! *  reach the next output time, do the output required, and continue taking
! *  integration steps and invoking output this way until the final time is
! *  reached, which triggers a `break' to jump out of the infinite loop set up
! *  with `while(true)'.
! *-----------------------------------------------------------------------------
! */
!
subroutine evolve(mass, pos, vel, &
             n,  t, dt_param, dt_dia, dt_out, &
            dt_tot,  init_out,  x_flag)
	use numz
	use dim
	use face, only:get_acc_jerk_pot_coll,write_diagnostics &
	               ,put_snapshot,evolve_step
	implicit none
    real(b8) mass(0:), pos(0:,0:)
    real(b8) vel(0:,0:)
    integer  n
    real(b8) t, dt_param, dt_dia, dt_out, dt_tot
    logical init_out,x_flag

    real(b8), allocatable:: acc(:,:),jerk(:,:)
    real(b8) epot                    ! potential energy of the n-body system
    real(b8) coll_time               ! collision (close encounter) time scale
    real(b8) einit                   ! initial total energy of the system
    integer nsteps                   ! number of integration time steps completed
    real(b8) dt
    real(b8) t_dia,t_out,t_end

    write(*,*) " Starting a Hermite integration for a " , n ,"-body system"
    write(*,*) " from time t = " , t
    write(*,*) " with time step control parameter dt_param = " , dt_param
    write(*,*) " until time " , t + dt_tot
    write(*,*) " with diagnostics output interval dt_dia = ", dt_dia
    write(*,*) " and snapshot output interval dt_out = ",dt_out ,"."

    allocate(acc(0:n-1,0:ndim-1))      ! accelerations and jerks
    allocate(jerk(0:n-1,0:ndim-1))     ! for all particles

    call get_acc_jerk_pot_coll(mass, pos, vel, acc, jerk, n, epot, coll_time)

    nsteps = 0               ! number of integration time steps completed

    call write_diagnostics(mass, pos, vel, acc, jerk, n, t, epot, &
                           nsteps, einit, .true., x_flag)
    if (init_out) then                 ! flag for initial output
        call put_snapshot(mass, pos, vel, n, t)
    endif
    t_dia = t + dt_dia           ! next time for diagnostics output
    t_out = t + dt_out           ! next time for snapshot output
    t_end = t + dt_tot           ! final time, to finish the integration

    do
        do while (t .lt. t_dia .and. t .lt. t_out .and. t .lt. t_end)
            dt = dt_param * coll_time
            call evolve_step(mass, pos, vel, acc, jerk, n, t, dt, epot, coll_time)
            nsteps=nsteps+1
        end do
        if (t .ge. t_dia)then
            call write_diagnostics(mass, pos, vel, acc, jerk, n, t, epot, nsteps, &
                                   einit, .false., x_flag)
            t_dia = t_dia + dt_dia
        endif
        if (t .ge. t_out)then
            call put_snapshot(mass, pos, vel, n, t)
            t_out = t_out + dt_out
        endif
        if (t .ge. t_end)exit
    end do

    deallocate(acc)
    deallocate(jerk)
end subroutine evolve

!/*-----------------------------------------------------------------------------
! *  put_snapshot  --  writes a single snapshot on the output stream cout.
! *
! *  note: unlike get_snapshot(), put_snapshot handles particle number and time
! *-----------------------------------------------------------------------------
! */
!
subroutine put_snapshot(mass, pos, vel, n, t)
 	use numz
	use dim
	implicit none
    real(b8) mass(0:), pos(0:,0:)
    real(b8) vel(0:,0:)
    integer  n
    real(b8) t
    integer i,k

                               ! form1 = full double precision

    write(*,*)n                         ! N, total particle number
    write(*,form1) t                                  ! current time
    do i=0,n-1
      write(*,form1,advance="no")mass(i)              ! mass of particle i
      do k=0,ndim-1
            write(*,form1,advance="no") pos(i,k)      ! position of particle i
      enddo
      do k=0,ndim-1
            write(*,form1,advance="no") vel(i,k)      ! velocity of particle i
      enddo
     write(*,*)
    enddo

end subroutine put_snapshot

!/*-----------------------------------------------------------------------------
! *  write_diagnostics  --  writes diagnostics on the error stream cerr:
! *                         current time; number of integration steps so far;
! *                         kinetic, potential, and total energy; absolute and
! *                         relative energy errors since the start of the run.
! *                         If x_flag (x for eXtra data) is true, all internal
! *                         data are dumped for each particle (mass, position,
! *                         velocity, acceleration, and jerk).
! *
! *  note: the kinetic energy is calculated here, while the potential energy is
! *        calculated in the function get_acc_jerk_pot_coll().
! *-----------------------------------------------------------------------------
! */
!
    subroutine write_diagnostics(mass, pos, &
                       vel, acc, &
                       jerk, n, t, epot, &
                       nsteps, einit,init_flag, &
                       x_flag)
	use numz
	use dim
	implicit none
    real(b8) mass(0:), pos(0:,0:)
    real(b8) vel(0:,0:),  acc(0:,0:)
    real(b8) jerk(0:,0:)
    integer  n
    real(b8) t
    real(b8) epot
    integer nsteps
    real(b8) einit
    logical init_flag,x_flag
    real(b8) ekin                        ! kinetic energy of the n-body system
    real(b8) etot
    integer i,k
    ekin = 0.0_b8
                           ! form1 =full double precision

    do i=0,n-1
       do k=0,ndim-1
            ekin = ekin+ 0.5 * mass(i) * vel(i,k) * vel(i,k)
       enddo
    enddo
    etot = ekin + epot             ! total energy of the n-body system

    if (init_flag) then                       ! at first pass, pass the initial
        einit = etot                    ! energy back to the calling function
    endif
    write(cerr,*)"at time t = " ,t, " , after " ,nsteps," steps"
    write(cerr,*)"E_kin = " ,ekin, "E_pot = " ,epot," , E_tot = ",etot
    write(cerr,*)"absolute energy error: E_tot - E_init = ",etot - einit
    write(cerr,*)"relative energy error: (E_tot - E_init) / E_init = ", &
                 (etot - einit) / einit

    if (x_flag)then
         write(cerr,*)"for debugging purposes, here is",&
                      " the internal data representation:\n"
        do i=0,n-1
            write(cerr,*)"   internal data for particle ", i+1 ," : "
            write(cerr,form1,advance="no")mass(i)
            do k=0,ndim-1
              write(*,form1,advance="no") pos(i,k)
            enddo
            do k=0,ndim-1
              write(*,form1,advance="no") vel(i,k)
            enddo
            do k=0,ndim-1
              write(*,form1,advance="no") acc(i,k)
            enddo
            do k=0,ndim-1
              write(*,form1,advance="no") jerk(i,k)
            enddo
            write(cerr,*)
        enddo
    endif
end subroutine write_diagnostics!/*-----------------------------------------------------------------------------
! *  read_options  --  reads the command line options, and implements them.
! *
! *  note: when the help option -h is invoked, the return value is set to false,
! *        to prevent further execution of the main program; similarly, if an
! *        unknown option is used, the return value is set to false.
! *-----------------------------------------------------------------------------
! */
subroutine read_options(flag, dt_param, dt_dia, dt_out, dt_tot, i_flag,x_flag)
	use numz
	use dim
	implicit none
    logical flag
    real(b8) dt_param,dt_dia,dt_out,dt_tot
    logical i_flag,x_flag
    namelist /options/ dt_param,dt_dia,dt_out,dt_tot,i_flag,x_flag
    integer istat,jstat
    flag=.true.
    open(unit=19,file="options.in",status="old",iostat=istat)
    if(istat .eq. 0) then
      read(19,nml=options,iostat=jstat)
      close(19)
      if(jstat .ne. 0)flag=.false.
      return
    endif
    write(cerr,nml=options)
end subroutine read_options

!/*-----------------------------------------------------------------------------
! *  get_snapshot  --  reads a single snapshot from the input stream cin.
! *
! *  note: in this implementation, only the particle data are read in, and it
! *        is left to the main program to first read particle number and time
! *-----------------------------------------------------------------------------
! */
!
subroutine get_snapshot(mass,  pos,  vel, n)
	use numz
	use dim
	implicit none
    real(b8) mass(0:), pos(0:,0:)
    real(b8) vel(0:,0:)
    integer n
    integer i,k
    write(*,*)lbound(mass),ubound(mass)
    do i=0,n-1
      read(*,*)mass(i),(pos(i,k),k=0,ndim-1),(vel(i,k),k=0,ndim-1)
    enddo
end subroutine get_snapshot!              ! Time-stamp: <2002-01-18 21:51:36 piet>
!             !================================================================
!            !                                                                |
!           !           /__----__                         ........            |
!          !       .           \                     ....:        :.          |
!         !       :                 _\|/_         ..:                         |
!        !       :                   /|\         :                     _\|/_  |
!       !  ___   ___                  _____                      ___    /|\   |
!      !  /     |   \    /\ \     / |   |   \  / |        /\    |   \         |
!     !  |   __ |___/   /  \ \   /  |   |    \/  |       /  \   |___/         |
!    !   |    | |  \   /____\ \ /   |   |    /   |      /____\  |   \     \/  |
!   !     \___| |   \ /      \ V    |   |   /    |____ /      \ |___/     |   |
!  !                                                                      /   |
! !              :                       _/|     :..                    |/    |
!!                :..               ____/           :....          ..         |
!/*   o   !          :.    _\|/_     /                   :........:           |
! *  O  `!\                 /|\                                               |
! *  |     /\                                                                  |
! *=============================================================================
! *
! *  nbody_sh1.C:  an N-body integrator with a shared but variable time step
! *                (the same for all particles but changing in time), using
! *                the Hermite integration scheme.
! *
! *                ref.: Hut, P., Makino, J. & McMillan, S., 1995,
! *                      Astrophysical Journal Letters 443, L93-L96.
! *
! *  note: in this first version, all functions are included in one file,
! *        without any use of a special library or header files.
! *_____________________________________________________________________________
! *
! *  usage: nbody_sh1 [-h (for help)] [-d step_size_control_parameter]
! *                   [-e diagnostics_interval] [-o output_interval]
! *                   [-t total_duration] [-i (start output at t = 0)]
! *                   [-x (extra debugging diagnostics)]
! *
! *         "step_size_control_parameter" is a coefficient determining the
! *            the size of the shared but variable time step for all particles
! *
! *         "diagnostics_interval" is the time between output of diagnostics,
! *            in the form of kinetic, potential, and total energy with the
! *            -x option, a dump of the internal particle data is made as well
! *
! *         "output_interval" is the time between successive snapshot outputs
! *
! *         "total_duration" is the integration time, until the program stops
! *
! *         Input and output are written from the standard i/o streams.  Since
! *         all options have sensible defaults, the simplest way to run the code
! *         is by only specifying the i/o files for the N-body snapshots:
! *
! *            nbody_sh1 < data.in > data.out
! *
! *         The diagnostics information will then appear on the screen.
! *         To capture the diagnostics information in a file, capture the
! *         standard error stream as follows:
! *
! *            (nbody_sh1 < data.in > data.out) >& data.err
! *
! *  Note: if any of the times specified in the -e, -o, or -t options are not an
! *        an integer multiple of "step", output will occur slightly later than
! *        predicted, after a full time step has been taken.  And even if they
! *        are integer multiples, round-off error may induce one extra step.
! *_____________________________________________________________________________
! *
! *  External data format:
! *
! *     The program expects input of a single snapshot of an N-body system,
! *     in the following format: the number of particles in the snapshot n
! *     the time t mass mi, position ri and velocity vi for each particle i,
! *     with position and velocity given through their three Cartesian
! *     coordinates, divided over separate lines as follows:
! *
! *                      n
! *                      t
! *                      m1 r1_x r1_y r1_z v1_x v1_y v1_z
! *                      m2 r2_x r2_y r2_z v2_x v2_y v2_z
! *                      ...
! *                      mn rn_x rn_y rn_z vn_x vn_y vn_z
! *
! *     Output of each snapshot is written according to the same format.
! *
! *  Internal data format:
! *
! *     The data for an N-body system is stored internally as a 1-dimensional
! *     array for the masses, and 2-dimensional arrays for the positions,
! *     velocities, accelerations and jerks of all particles.
! *_____________________________________________________________________________
! *
! *    version 1:  Jan 2002   Piet Hut, Jun Makino
! *_____________________________________________________________________________
! */
!

!/*-----------------------------------------------------------------------------
! *  main  --  reads option values, reads a snapshot, and launches the
! *            integrator
! *-----------------------------------------------------------------------------
! */

program ell
	use numz
	use dim
	use face, only:read_options,get_snapshot,evolve
	implicit none
    logical flag

    real(b8)  dt_param        ! control parameter to determine time step size
    real(b8)  dt_dia          ! time interval between diagnostics output
    real(b8)  dt_out          ! time interval between output of snapshots
    real(b8)  dt_tot          ! duration of the integration
    logical  init_out         ! if true: snapshot output with start at t = 0
                              !          with an echo of the input snapshot
    logical  x_flag           ! if true: extra debugging diagnostics output
    logical  read_flag
    integer n
    real(b8) t
    real(b8),allocatable :: mass(:),pos(:,:),vel(:,:)

	dt_param = 0.03
	dt_dia = 1
	dt_out = 1
	dt_tot = 10
	init_out = .false.

	x_flag = .false.

    call read_options(read_flag, dt_param, dt_dia, dt_out, &
                      dt_tot, init_out, x_flag)
                             ! halt criterion detected by read_options()
    if (read_flag .eq. .false.)stop

    read(*,*)n             ! N, number of particles in the N-body system

    read(*,*)t             ! time

    allocate(mass(0:n-1))             ! masses for all particles
    allocate(pos(0:n-1,0:NDIM-1))     ! positions for all particles
    allocate(vel(0:n-1,0:NDIM-1))     ! velocities for all particles

    call get_snapshot(mass, pos, vel, n)

    call evolve(mass, pos, vel, n, t, dt_param, dt_dia, &
                dt_out, dt_tot, init_out, x_flag)

    deallocate(mass)
    deallocate(pos)
    deallocate(vel)
end program

max_stars_repo_name: timkphd/examples

max_issues_repo_name: timkphd/examples

max_forks_repo_name: timkphd/examples


text:
      SUBROUTINE PLTOUT
     &           (NSD,   NNW_SD, NW_IPAR,
     &            CRD,   PHI,    PHN,
     &            CRD_T, PHI_T,  PHN_T)
C ---------------------------------------------------------------------------
C
C ---------------------------------------------------------------------------
      IMPLICIT NONE
      INCLUDE 'knd_params.inc'
      INCLUDE 'usr_params.inc'
      INCLUDE 'fle_params.inc'
      INCLUDE 'net_params.inc'
      INCLUDE 'net_funcs.inc'
      INCLUDE 'tme_funcs.inc'
C
      INTEGER(KIND=IK) NSD
      INTEGER(KIND=IK) NNW_SD (*)
      INTEGER(KIND=IK) NW_IPAR (N_IP, *)
      REAL   (KIND=RK) CRD (4, *)
      REAL   (KIND=RK) PHI (2, *)
      REAL   (KIND=RK) PHN (2, *)
      REAL   (KIND=RK) CRD_T (4, *)
      REAL   (KIND=RK) PHI_T (2, *)
      REAL   (KIND=RK) PHN_T (2, *)
C
      INTEGER(KIND=IK) ISD, INW, JNW, NNW, IGP, NGP, PLT, BCT, CDF
      REAL   (KIND=RK) T, PERIOD
      REAL   (KIND=RK) NRM (6)
C
      PERIOD = GET_PERIOD ()
      T      = GET_TIME   ()
      WRITE (PLT_CRD, '(''# T = '',F8.4, '' * PERIOD'')') T / PERIOD
      WRITE (PLT_PHI, '(''# T = '',F8.4, '' * PERIOD'')') T / PERIOD
      WRITE (PLT_PHN, '(''# T = '',F8.4, '' * PERIOD'')') T / PERIOD
      WRITE (PLT_NRM, 101) 'T_REL',
     &                     'PHI_MX', 'PHI_L1', 'PHI_L2',
     &                     'PHN_MX', 'PHN_L1', 'PHN_L2'
C
      INW = 1
      IGP = 1
C
      DO ISD = 1, NSD
         NNW = NNW_SD (ISD)
         DO JNW = INW, INW + NNW - 1
            NGP = GET_NGP (NW_IPAR (1, JNW))
            BCT = GET_BCT (NW_IPAR (1, JNW))
            PLT = GET_PLT (NW_IPAR (1, JNW))
            IF (PLT .NE. 0) THEN
               CALL PLT_NW
     &              (NGP, ABS (BCT), PLT,
     &               T,
     &               CRD   (1, IGP), PHI   (1, IGP), PHN   (1, IGP),
     &               CRD_T (1, IGP), PHI_T (1, IGP), PHN_T (1, IGP),
     &               NRM)
               WRITE (PLT_NRM, 121) T / PERIOD, NRM
            END IF
            IGP = IGP + NGP
         END DO
         INW = INW + NNW
      END DO
      WRITE (PLT_NRM, *)
C
      RETURN
  101 FORMAT (1X, '#', A13, 6A12)
  121 FORMAT (1X, E14.6, 10E12.4)
      END
      SUBROUTINE PLT_NW
     &           (NGP,   BCT,   PLT,
     &            T,
     &            CRD,   PHI,   PHN,
     &            CRD_T, PHI_T, PHN_T, NRM)
C ---------------------------------------------------------------------------
C
C ---------------------------------------------------------------------------
      IMPLICIT NONE
      INCLUDE 'knd_params.inc'
      INCLUDE 'usr_params.inc'
      INCLUDE 'anl_params.inc'
      INCLUDE 'net_params.inc'
      INCLUDE 'bct_params.inc'
      INCLUDE 'fle_params.inc'
C
      INTEGER(KIND=IK) NGP
      INTEGER(KIND=IK) BCT
      INTEGER(KIND=IK) PLT
      REAL   (KIND=RK) T
      REAL   (KIND=RK) CRD (4, NGP)
      REAL   (KIND=RK) PHI (2, *)
      REAL   (KIND=RK) PHN (2, *)
      REAL   (KIND=RK) CRD_T (4, *)
      REAL   (KIND=RK) PHI_T (2, *)
      REAL   (KIND=RK) PHN_T (2, *)
      REAL   (KIND=RK) NRM   (3, *)
C
      INTEGER(KIND=IK) I, J
      REAL   (KIND=RK) X_XI, Z_XI, J_I, PROD, DIFF
      REAL   (KIND=RK) W (NGP, 8)
      REAL   (KIND=RK) U_L (2), U_G (2)
C
      INTEGER(KIND=IK) IDAMAX
      EXTERNAL         IDAMAX
      REAL   (KIND=RK) DNRM2, DASUM
      EXTERNAL         DNRM2, DASUM
C
      IF (BCT .EQ. BCT_LINBERNOU) THEN
         DO I = 1, NGP
            W   (I, 8) = CRD (2, I)
            CRD (2, I) = 0.0D0
         END DO
      END IF
      CALL ANALYTIC (NGP, T, ANL_PAR, CRD, W)
      IF (BCT .EQ. BCT_LINBERNOU) THEN
         DO I = 1, NGP
            CRD (2, I) = W (I, 8)
         END DO
      END IF
C
       DO I = 1, NGP
         X_XI  = CRD (3, I)
         Z_XI  = CRD (4, I)
         J_I   = 1.0D0 / SQRT (X_XI ** 2 + Z_XI ** 2)
         PROD  = 2.0D0 * X_XI * Z_XI
         DIFF  = (X_XI + Z_XI) * (X_XI - Z_XI)
         W (I, 1) = PHI (1, I) -  W (I,1)
         W (I, 2) = PHI (2, I) - (W (I,3) * X_XI + W (I,4) * Z_XI)
         W (I, 3) = PHN (1, I) - (W (I,4) * X_XI - W (I,3) * Z_XI) * J_I
         W (I, 4) = PHN (2, I) - (W (I,6) * DIFF - W (I,5) * PROD) * J_I
         W (I, 5) = CRD (2, I) -  W (I,7)
         W (I, 6) = ABS (W (I, 1))
         W (I, 7) = ABS (W (I, 3))
         U_L (1) = PHI (2, I) * J_I
         U_L (2) = PHN (1, I)
         U_G (1) = J_I * (CRD (3, I) * U_L (1) - CRD (4, I) * U_L (2))
         U_G (2) = J_I * (CRD (3, I) * U_L (2) + CRD (4, I) * U_L (1))
         WRITE (PLT_CRD, 1001) (CRD (J,I), J = 1, 4), U_G
         WRITE (PLT_PHI, 1001)  PHI (1,I), PHI (2,I), W (I,1), W (I,2)
         WRITE (PLT_PHN, 1001)  PHN (1,I), PHN (2,I), W (I,3), W (I,4)
      END DO
      WRITE (PLT_CRD, *)
      WRITE (PLT_PHI, *)
      WRITE (PLT_PHN, *)
C
      NRM (1, 1) = W (IDAMAX (NGP, W (1, 6), 1), 6)
      NRM (2, 1) =     DASUM (NGP, W (1, 6), 1)
      NRM (3, 1) =     DNRM2 (NGP, W (1, 6), 1)
      NRM (1, 2) = W (IDAMAX (NGP, W (1, 7), 1), 7)
      NRM (2, 2) =     DASUM (NGP, W (1, 7), 1)
      NRM (3, 2) =     DNRM2 (NGP, W (1, 7), 1)
C
      RETURN
 1001 FORMAT (1X,  10E16.8)
      END

max_stars_repo_name: meggermo/waterwaves

max_issues_repo_name: meggermo/waterwaves

max_forks_repo_name: meggermo/waterwaves


text:
c
c
c     ##############################################################
c     ##  COPYRIGHT (C) 1997 by Rohit Pappu & Jay William Ponder  ##
c     ##                   All Rights Reserved                    ##
c     ##############################################################
c
c     ##################################################################
c     ##                                                              ##
c     ##  program pssrot  --  torsional potential smoothing & search  ##
c     ##                                                              ##
c     ##################################################################
c
c
c     "pssrot" implements the potential smoothing plus search method
c     for global optimization in torsional space
c
c     literature reference:
c
c     J. Kostrowicki and H. A. Scheraga, "Application of the Diffusion
c     Equation Method for Global Optimization to Oligopeptides", Journal
c     of Physical Chemistry, 96, 7442-7449 (1992)
c
c     S. Nakamura, H. Hirose, M. Ikeguchi and J. Doi, "Conformational
c     Energy Minimization Using a Two-Stage Method", Journal of Physical
c     Chemistry, 99, 8374-8378 (1995)
c
c
      program pssrot
      implicit none
      include 'sizes.i'
      include 'atoms.i'
      include 'files.i'
      include 'inform.i'
      include 'iounit.i'
      include 'math.i'
      include 'omega.i'
      include 'refer.i'
      include 'warp.i'
      include 'zcoord.i'
      integer i,k
      integer ixyz,next
      integer npoint,neigen
      integer lext,freeunit
      real*8 minimum,grdmin
      real*8 pssrot1,rms
      real*8 srchmax
      real*8 deform0,ratio
      real*8, allocatable :: xx(:)
      logical exist
      logical use_local
      character*1 answer
      character*7 ext
      character*120 xyzfile
      character*120 record
      character*120 string
      external pssrot1
      external optsave
c
c
c     set up the structure, mechanics calculation and smoothing
c
      call initial
      call getint
      use_smooth = .true.
      use_dem = .true.
      call mechanic
      call initrot
c
c     convert to Cartesian coordinates and save the initial set
c
      call makexyz
      call makeref (1)
c
c     set maximum deformation value and disable coordinate dumps
c
      deform0 = deform
      iwrite = 0
c
c     get the number of points along the deformation schedule
c
      npoint = -1
      call nextarg (string,exist)
      if (exist)  read (string,*,err=10,end=10)  npoint
   10 continue
      if (npoint .lt. 0) then
         write (iout,20)
   20    format (/,' Enter the Number of Steps for the PSS Schedule',
     &              ' [100] :  ',$)
         read (input,30)  npoint
   30    format (i10)
         if (npoint .le. 0)  npoint = 100
      end if
c
c     decide whether to use the local search procedure
c
      use_local = .false.
      call nextarg (answer,exist)
      if (.not. exist) then
         write (iout,40)
   40    format (/,' Use Local Search to Explore the Smoothing Levels',
     &              ' [N] :  ',$)
         read (input,50)  record
   50    format (a120)
         next = 1
         call gettext (record,answer,next)
      end if
      call upcase (answer)
      if (answer .eq. 'Y')  use_local = .true.
c
c     get the number of eigenvectors to use for the local search
c
      if (use_local) then
         neigen = -1
         call nextarg (string,exist)
         if (exist)  read (string,*,err=60,end=60)  neigen
   60    continue
         if (neigen .le. 0) then
            write (iout,70)
   70       format (/,' Enter the Number of Directions for Local',
     &                 ' Search [5] :  ',$)
            read (input,80)  neigen
   80       format (i10)
            if (neigen .le. 0)  neigen = 5
         end if
      end if
c
c     get the maximal smoothing level for use of local search
c
      if (use_local) then
         srchmax = -1.0d0
         call nextarg (string,exist)
         if (exist)  read (string,*,err=90,end=90)  srchmax
   90    continue
         if (srchmax .lt. 0.0d0) then
            write (iout,100)
  100       format (/,' Enter the Largest Smoothing Value for Local',
     &                 ' Search [5.0] :  ',$)
            read (input,110)  srchmax
  110       format (f20.0)
            if (srchmax .lt. 0.0d0)  srchmax = 5.0d0
         end if
      end if
c
c     get the termination criterion as RMS gradient per atom
c
      grdmin = -1.0d0
      call nextarg (string,exist)
      if (exist)  read (string,*,err=120,end=120)  grdmin
  120 continue
      if (grdmin .le. 0.0d0) then
         write (iout,130)
  130    format (/,' Enter RMS Gradient per Atom Criterion',
     &              ' [0.0001] :  ',$)
         read (input,140)  grdmin
  140    format (f20.0)
      end if
      if (grdmin .le. 0.0d0)  grdmin = 0.0001d0
c
c     perform dynamic allocation of some local arrays
c
      allocate (xx(nomega))
c
c     perform PSS iteration by looping over smoothed surfaces
c
      do k = 0, 2*npoint
         ratio = 1.0d0 - dble(abs(npoint-k))/dble(npoint)
         deform = deform0 * ratio**3
c
c     translate the initial coordinates
c
         do i = 1, nomega
            xx(i) = dihed(i)
         end do
c
c     make the call to the variable metric optimization routine
c
         iprint = 1
         call ocvm (nomega,xx,minimum,grdmin,pssrot1,optsave)
c
c     untranslate the final coordinates for each atom
c
         do i = 1, nomega
            dihed(i) = xx(i)
            ztors(zline(i)) = dihed(i) * radian
         end do
c
c     use normal mode local search to explore adjacent minima
c
         if (use_local) then
            if (k.ge.npoint .and. deform.le.srchmax)
     &         call moderot (neigen,minimum,grdmin)
         end if
c
c     write out final energy function value and smoothing level
c
         write (iout,150)  minimum,deform
  150    format (/,' Final Function Value and Deformation :',2f15.4)
c
c     get Cartesian coordinates and superimpose on reference
c
         call makexyz
         call impose (n,xref,yref,zref,n,x,y,z,rms)
c
c     write the coordinates of the current minimum to a file
c
         lext = 3
         call numeral (k,ext,lext)
         ixyz = freeunit ()
         xyzfile = filename(1:leng)//'.'//ext(1:lext)
         call version (xyzfile,'new')
         open (unit=ixyz,file=xyzfile,status='new')
         call prtxyz (ixyz)
         close (unit=ixyz)
      end do
c
c     perform deallocation of some local arrays
c
      deallocate (xx)
c
c     perform any final tasks before program exit
c
      call final
      end
c
c
c     ################################################################
c     ##                                                            ##
c     ##  function pssrot1  --  energy and gradient values for PSS  ##
c     ##                                                            ##
c     ################################################################
c
c
c     "pssrot1" is a service routine that computes the energy and
c     gradient during PSS global optimization in torsional space
c
c
      function pssrot1 (xx,g)
      implicit none
      include 'sizes.i'
      include 'math.i'
      include 'omega.i'
      include 'zcoord.i'
      integer i
      real*8 pssrot1,e
      real*8 xx(*)
      real*8 g(*)
      real*8, allocatable :: derivs(:)
c
c
c     translate optimization variables into dihedrals
c
      do i = 1, nomega
         dihed(i) = xx(i)
         ztors(zline(i)) = dihed(i) * radian
      end do
c
c     perform dynamic allocation of some local arrays
c
      allocate (derivs(nomega))
c
c     compute and store the energy and gradient
c
      call makexyz
      call gradrot (e,derivs)
      pssrot1 = e
c
c     store torsional gradient as optimization gradient
c
      do i = 1, nomega
         g(i) = derivs(i)
      end do
c
c     perform deallocation of some local arrays
c
      deallocate (derivs)
      return
      end
c
c
c     ##############################################################
c     ##                                                          ##
c     ##  subroutine moderot  --  torsional local search for PSS  ##
c     ##                                                          ##
c     ##############################################################
c
c
      subroutine moderot (neigen,minimum,grdmin)
      implicit none
      include 'sizes.i'
      include 'iounit.i'
      include 'math.i'
      include 'omega.i'
      include 'zcoord.i'
      integer i,k,neigen
      integer ndoi,nsearch
      real*8 minimum,grdmin
      real*8 eps,minref,minbest
      real*8, allocatable :: step(:)
      real*8, allocatable :: eigen(:)
      real*8, allocatable :: zorig(:)
      real*8, allocatable :: zbest(:)
      real*8, allocatable :: vects(:,:)
      logical done
c
c
c     perform dynamic allocation of some local arrays
c
      allocate (step(nomega))
      allocate (eigen(nomega))
      allocate (zorig(nomega))
      allocate (zbest(nomega))
      allocate (vects(nomega,nomega))
c
c     set parameters related to the local search procedure
c
      done = .false.
      eps = 0.0001d0
      minref = minimum
      minbest = minimum
      ndoi = 0
      do k = 1, nomega
         zorig(k) = ztors(zline(k))
      end do
c
c     find local minimum along each of the steepest directions
c
      do while (.not. done)
         ndoi = ndoi + 1
         write (iout,10)  ndoi,minref
   10    format (/,' Torsional Mode Search :',5x,'Iteration',i4,
     &              6x,'Energy',f12.4,/)
         call makexyz
         call eigenrot (eigen,vects)
c
c     search both directions along each eigenvector in turn
c
         nsearch = 0
         do i = 1, neigen
            do k = 1, nomega
               step(k) = vects(k,nomega-i+1)
               ztors(zline(k)) = zorig(k)
            end do
            nsearch = nsearch + 1
            call climbrot (nsearch,minimum,step,grdmin)
            if (minimum .lt. minbest) then
               minbest = minimum
               do k = 1, nomega
                  zbest(k) = ztors(zline(k))
               end do
            end if
            do k = 1, nomega
               step(k) = -vects(k,nomega-i+1)
               ztors(zline(k)) = zorig(k)
            end do
            nsearch = nsearch + 1
            call climbrot (nsearch,minimum,step,grdmin)
            if (minimum .lt. minbest) then
               minbest = minimum
               do k = 1, nomega
                  zbest(k) = ztors(zline(k))
               end do
            end if
         end do
c
c     check for convergence of the local search procedure
c
         if (minbest .lt. minref-eps) then
            done = .false.
            minref = minbest
            do k = 1, nomega
               zorig(k) = zbest(k)
            end do
         else
            done = .true.
            minimum = minref
            do k = 1, nomega
               dihed(k) = zorig(k) / radian
               ztors(zline(k)) = zorig(k)
            end do
         end if
      end do
c
c     perform deallocation of some local arrays
c
      deallocate (step)
      deallocate (eigen)
      deallocate (zorig)
      deallocate (zbest)
      deallocate (vects)
      return
      end
c
c
c     ###############################################################
c     ##                                                           ##
c     ##  subroutine eigenrot  --  torsional Hessian eigenvectors  ##
c     ##                                                           ##
c     ###############################################################
c
c
      subroutine eigenrot (eigen,vects)
      implicit none
      include 'sizes.i'
      include 'atoms.i'
      include 'omega.i'
      integer i,j,ihess
      real*8 eigen(*)
      real*8, allocatable :: matrix(:)
      real*8 vects(nomega,*)
      real*8, allocatable :: hrot(:,:)
c
c
c     perform dynamic allocation of some local arrays
c
      allocate (matrix(nomega*(nomega+1)/2))
      allocate (hrot(nomega,nomega))
c
c     compute the Hessian in torsional space
c
      call hessrot ('FULL',hrot)
c
c     place Hessian elements into triangular form
c
      ihess = 0
      do i = 1, nomega
         do j = i, nomega
            ihess = ihess + 1
            matrix(ihess) = hrot(i,j)
         end do
      end do
c
c     diagonalize the Hessian to obtain eigenvalues
c
      call diagq (nomega,nomega,matrix,eigen,vects)
c
c     perform deallocation of some local arrays
c
      deallocate (matrix)
      deallocate (hrot)
      return
      end
c
c
c     ################################################################
c     ##                                                            ##
c     ##  subroutine climbrot  --  minimum from a PSS local search  ##
c     ##                                                            ##
c     ################################################################
c
c
      subroutine climbrot (nsearch,minimum,step,grdmin)
      implicit none
      include 'sizes.i'
      include 'iounit.i'
      include 'math.i'
      include 'omega.i'
      include 'zcoord.i'
      integer maxstep
      parameter (maxstep=500)
      integer i,nsearch
      integer kstep,nstep
      real*8 minimum,grdmin
      real*8 energy
      real*8 big,small,size
      real*8 estep(0:maxstep)
      real*8 step(*)
      logical done
c
c
c     set the maximum number of steps and the step size
c
      done = .false.
      big = 1.0d10
      small = -1.0d5
      minimum = big
      kstep = 0
      nstep = 65
      size = 0.1d0 * radian
      do i = 1, nomega
         step(i) = size * step(i)
      end do
c
c     scan the search direction for a minimization candidate
c
      do while (.not. done)
         if (kstep .ne. 0) then
            do i = 1, nomega
               ztors(zline(i)) = ztors(zline(i)) + step(i)
            end do
         end if
         call makexyz
         estep(kstep) = energy ()
         if (kstep .ge. 2) then
            if (estep(kstep) .lt. estep(kstep-2) .and.
     &          estep(kstep-1) .lt. estep(kstep-2)) then
               done = .true.
               do i = 1, nomega
                  ztors(zline(i)) = ztors(zline(i)) - step(i)
               end do
               call makexyz
               call localrot (minimum,grdmin)
               if (minimum .ge. small) then
                  write (iout,10)  nsearch,kstep-1,minimum
   10             format (4x,'Search Direction',i4,10x,'Step',
     &                       i6,10x,f12.4)
               else
                  minimum = big
                  write (iout,20)  nsearch
   20             format (4x,'Search Direction',i4,36x,'------')
               end if
            end if
         end if
         if (kstep.ge.nstep .and. .not.done) then
            done = .true.
            write (iout,30)  nsearch
   30       format (4x,'Search Direction',i4,36x,'------')
         end if
         kstep = kstep + 1
      end do
      return
      end
c
c
c     ##############################################################
c     ##                                                          ##
c     ##  subroutine localrot  --  PSS local search optimization  ##
c     ##                                                          ##
c     ##############################################################
c
c
c     "localrot" is used during the PSS local search procedure
c     to perform a torsional space energy minimization
c
c
      subroutine localrot (minimum,grdmin)
      implicit none
      include 'sizes.i'
      include 'inform.i'
      include 'minima.i'
      include 'math.i'
      include 'omega.i'
      include 'zcoord.i'
      integer i,oldprt
      real*8 minimum
      real*8 grdmin
      real*8 pssrot1
      real*8, allocatable :: xx(:)
      logical oldverb
      external pssrot1
      external optsave
c
c
c     perform dynamic allocation of some local arrays
c
      allocate (xx(nomega))
c
c     translate the coordinates of each atom
c
      do i = 1, nomega
         dihed(i) = ztors(zline(i)) / radian
         xx(i) = dihed(i)
      end do
c
c     make the call to the optimization routine
c
      oldverb = verbose
      oldprt = iprint
      verbose = .false.
      iprint = 0
      call ocvm (nomega,xx,minimum,grdmin,pssrot1,optsave)
      verbose = oldverb
      iprint = oldprt
c
c     untranslate the final coordinates for each atom
c
      do i = 1, nomega
         dihed(i) = xx(i)
         ztors(zline(i)) = dihed(i) * radian
      end do
c
c     perform deallocation of some local arrays
c
      deallocate (xx)
      return
      end

max_stars_repo_name: htna/HCsbLib

max_issues_repo_name: htna/HCsbLib

max_forks_repo_name: htna/HCsbLib