text:
C     Concrete Stress-Strain Model Including Unloading/Reloading Model
C     Proposed by Sakai and Kawashima (Envelope Curve: Mander Model)
C                                      Coded by J. Sakai    2002. 2.14
C ====================================================================
C ====================================================================
C ====================================================================
      SUBROUTINE NLU014( NTABL, MTB, NTB, D, MNLI, MNLR, NP, P,
     1                   NHST, HST, NIHST, IHST,
     1                   DEPS, DSIG, TK, DEPSV,
     1                   NEL, IPRT, INIT, LDBG, IT, ILOOP)
      IMPLICIT REAL*8(A-H,O-Z)
cDEC$ ATTRIBUTES DLLEXPORT:: NLU014, /MLSVAL/
      COMMON /MLSVAL/ VAL1, VAL2, VALS
      DIMENSION NTB(NTABL),D(MTB,NTABL),NP(MNLI),P(MNLR),
     1          HST(NHST),IHST(NIHST)
c  ---------------------------------------------------------------------
c     Material Type: Confined Concrete
c     Envelope Curve: Mander Model
c     Unloading and Reloading Curve: Sakai-Kawashima Model
c
c     Material Parameters from MATNL1 Data
c        P(1) : Initial Slope or Initial Modulus of Elasticity       YMc
c        P(2) : Peak Stress                                          Sigcc
c        P(3) : Strain at Peak Stress                                EPScc
c
c     State Variables
c        HST( 1) : Initial Stiffness
c        HST( 2) : Strain
c        HST( 3) : Stress
c        HST( 4) : Delta Epsilon of the Previous Step                  DE0
c        HST( 5) : Stress at the Unloading Point on Envelope           Sigule
c        HST( 6) : Strain at the Unloading Point on Envelope           EPSule
c        HST( 7) : Stress at the Unloading Point                       Sigul
c        HST( 8) : Strain at the Unloading Point                       EPSul
c        HST( 9) : Plastic Strain                                      EPSpl
c        HST(10) : Stress at the Unloading Point after Un/Reloading    Suln
c        HST(11) : Slope of Reloding Branch                            YMrl
c        HST(12) : Tangential Slope of the Previous Step               YMtan
c        HST(13) : Stress at the Reloading Point                       Sigrl
c        HST(14) : Strain at the Reloading Point                       EPSrl
c        HST(15) : HST( 9) EPSpl of Previous Unloading                 EPSpl0
c        HST(16) : HST(10) Suln of Previous Reloading                  Suln0
c        HST(17) : Partial Reloading Ratio GammaRL                     GamRL
c
c        IHST(1) : Index of Current Status                      Jcon
c        Jcon= 1 : On the Envelope Curve
c        Jcon= 2 : Zero Stress
c        Jcon= 3 : On the Unloading Curve (Parabolic)
c        Jcon= 4 : On the Reloading Line (Linear)   0.2< Normalized Strain <1.0
c        Jcon= 5 : On the Reloading Curve (Parabolic) 0< Normalized Strain <0.2
c        Jcon= 6 : On the Reloading Curve (Parabolic) 0.75< Normalized Stress <1
c        Jcon= 7 : On the Reloading Line (Linear)     Normalized Strain >1.0
c        Jcon= 8 : Zero Stress (No Compressive Loading)
c        IHST(2) : Number of Cyclcles                           Ncyc
c        IHST(3) : Index of Previous Status                     Jcon0
c        IHST(4) : Number of Cyclcles of Previous Status        Ncyc0
c
c        DEPS   (IN)   : Incremental Strain
c        DSIG   (OUT)  : Incremental Stress
c        TK     (OUT)  : Current Status
c  =====================================================================
c  ---------------------------------------------------------------------
c   Stress Tensor at the Start of This Increment
      Sig0=HST(3)
c  ---------------------------------------------------------------------
c   Reference Variables
      EPS0=HST(2)
      DE=DEPS
c  ---------------------------------------------------------------------
c   Strain at the End of This Increment
      EPS1=EPS0+DE
c  ---------------------------------------------------------------------
c   User Defined Material Parameters
      YMc=P(1)
      Sigcc=P(2)
      EPScc=P(3)
c
      r=YMc*EPScc/(YMc*EPScc-Sigcc)
c  ---------------------------------------------------------------------
c   Initial State Variables (INIT=1)
      IF(INIT.EQ.1) THEN
      NHST=17
      NIHST=4
      HST(1)=TK
c
      DO 10 I=4,NHST
      HST(I)=0.0
   10 CONTINUE
c
      IHST(1)=8
      IHST(2)=1
      IHST(3)=8
      IHST(4)=1
c
c  ---------------------------------------------------------------------
c   State Variables at the Start of This Increment
      ELSE
      DE0=HST(4)
      Sigule=HST(5)
      EPSule=HST(6)
      Sigul=HST(7)
      EPSul=HST(8)
      EPSpl=HST(9)
      Suln=HST(10)
      YMrl=HST(11)
      YMtan0=HST(12)
      Sigrl=HST(13)
      EPSrl=HST(14)
      EPSpl0=HST(15)
      Suln0=HST(16)
      GamRL=HST(17)
c
      Jcon=IHST(1)
      Ncyc=IHST(2)
      Jcon0=IHST(3)
      Ncyc0=IHST(4)
c
      r=YMc*EPScc/(YMc*EPScc-Sigcc)
C ********************************************************************
      IF(DE.EQ.0.0) THEN
      Sig1=Sig0
      YMtan=YMtan0
      ELSE IF(DE.LT.0.0) THEN
c   Compressive Loading
      CALL COMPR14 (EPS1,Sig1,YMtan,Jcon,EPS0,Sig0,YMc,EPScc,Sigcc,r,
     1              EPSule,Sigule,EPSpl,Suln,YMrl,EPSul,Sigul,
     1              EPSrl,Sigrl,EPSpl0,Suln0,GamRL,
     1              Ncyc,Jcon0,Ncyc0)
      ELSE
c   Unloading (Tensile) Loading
      CALL TENSI14 (EPS1,Sig1,YMtan,Jcon,EPS0,Sig0,
     1              EPSule,Sigule,EPSpl,Suln,EPSul,Sigul,
     1              EPSrl,Sigrl,EPSpl0,Suln0,GamRL,
     1              Ncyc,Jcon0,Ncyc0)
      END IF
C ********************************************************************
c   New State Variables
      HST(2)=EPS1
      HST(3)=Sig1
      DSIG=Sig1-Sig0
      TK=YMtan
      HST(4)=DE
      HST(5)=Sigule
      HST(6)=EPSule
      HST(7)=Sigul
      HST(8)=EPSul
      HST(9)=EPSpl
      HST(10)=Suln
      HST(11)=YMrl
      HST(12)=YMtan
      HST(13)=Sigrl
      HST(14)=EPSrl
      HST(15)=EPSpl0
      HST(16)=Suln0
      HST(17)=GamRL
c
      IHST(1)=Jcon
      IHST(2)=Ncyc
      IHST(3)=Jcon0
      IHST(4)=Ncyc0
c  ---------------------------------------------------------------------
      END IF
      RETURN
      END
c
C ====================================================================
C ------------------------ SUBROUTINE PROGRAM ------------------------
C ------------------- Steel Stress-Strain Relation   -----------------
C ========================== COMPRESSION =============================
C ====================================================================
C
      SUBROUTINE COMPR14 (EPS1,Sig1,YMtan,Jcon,EPS0,Sig0,YMc,EPScc,
     1                    Sigcc,r,EPSule,Sigule,EPSpl,Suln,YMrl,EPSul,
     1                    Sigul,EPSrl,Sigrl,EPSpl0,Suln0,GamRL,
     1                    Ncyc,Jcon0,Ncyc0)
      IMPLICIT REAL*8(A-H,O-Z)
c
C *** ENVELOPE CURVE (Jcon=1) ****************************************
      IF(Jcon.EQ.1) THEN
      GO TO 20
C
C *** ZERO STRESS (Jcon=2) *******************************************
      ELSE IF(Jcon.EQ.2) THEN
        IF(EPS1.LT.EPSpl) THEN
        Jcon0=2
C
        YMrl=(Suln-0.1*Suln0)/(EPSule-(0.2*EPSule+0.8*EPSpl))
        Sigrl=0.0
        EPSrl=EPSpl
          IF(EPS1.LT.EPSule) THEN
          Jcon=7
          GO TO 10
          ELSE IF(EPS1.LE.0.2*EPSule+0.8*EPSpl) THEN
          Jcon=4
          Sig1=YMrl*(EPS1-EPSule)+Suln
          YMtan=YMrl
          ELSE
          Jcon=5
          Sig1=2.5*Suln0*((EPS1-EPSpl)/(EPSule-EPSpl))**2.0
          YMtan=5.0*Suln0*(EPS1-EPSpl)/((EPSule-EPSpl))**2.0
          END IF
        ELSE
C
C  When hysteresis does not reach the plastic strain
        Sig1=0.0
        YMtan=0.0
        END IF
C
C *** RELOADING FROM UNLOADING PATH (Jcon=3) *************************
      ELSE IF(Jcon.EQ.3) THEN
C
C  When partial unloading ratio is smaller than 0.25
        IF(ABS(Sig0-Sigul).LT.ABS(Sigule*0.25)) THEN
          IF(EPS1.LE.EPSul) THEN
C
C  When strain has already exceeded the previous unloaded strain
            IF(Jcon0.EQ.7) THEN
            Jcon=7
            Jcon0=3
              IF(Sigrl.EQ.0.0) THEN
              IF(GamRL.GE.0.7) THEN
              EPSpl=EPSpl0
              Suln=Suln0
              Ncyc=Ncyc0
              END IF
              END IF
            GO TO 10
            ELSE IF(Jcon0.EQ.4) THEN
            Jcon=4
            Jcon0=3
              IF(Sigrl.EQ.0.0) THEN
              IF(GamRL.GE.0.7) THEN
              EPSpl=EPSpl0
              Suln=Suln0
              Ncyc=Ncyc0
              END IF
              END IF
            GO TO 5
            ELSE
            Jcon=1
            Jcon0=3
            GO TO 20
            END IF
          ELSE
C
C  Otherwise, reloading along the previous unloading path
C  Jcon0 in not renewed here.
          Jcon=6
          Sig1=Sigul*((EPS1-EPSpl)/(EPSul-EPSpl))**2.0
          YMtan=2*Sigul*(EPS1-EPSpl)/((EPSul-EPSpl)**2.0)
          GO TO 30
          END IF
C
C  When partial unloading ratio is larger than 0.25
        ELSE
        Sigrl=Sig0
        EPSrl=EPS0
        YMrl=(Suln-Sigrl)/(EPSule-EPSrl)
        Jcon=4
        Jcon0=3
        GO TO 5
        END IF
C
C *** ON RELOADING PATH (Jcon=4) *************************************
      ELSE IF(Jcon.EQ.4) THEN
      GO TO 5
C
C *** ON RELOADING PATH, 0< Normalized Strain <0.2 (Jcon=5) **********
      ELSE IF(Jcon.EQ.5) THEN
        IF(EPS1.LT.EPSule) THEN
        Jcon=7
        GO TO 10
        ELSE IF(EPS1.LE.0.2*EPSule+0.8*EPSpl) THEN
        Jcon=4
        Sig1=YMrl*(EPS1-EPSule)+Suln
        YMtan=YMrl
        ELSE
        Sig1=2.5*Suln0*((EPS1-EPSpl)/(EPSule-EPSpl))**2.0
        YMtan=5.0*Suln0*(EPS1-EPSpl)/((EPSule-EPSpl)**2.0)
        END IF
C
C *** RELOADING ALONG THE PREVIOUS UNLOADING PATH (Jcon=6) ***********
      ELSE IF(Jcon.EQ.6) THEN
C
C  When strain has already exceeded the previous unloaded strain
        IF(EPS1.LE.EPSul) THEN
          IF(Jcon0.EQ.7) THEN
          Jcon=7
          Jcon0=3
            IF(Sigrl.EQ.0.0) THEN
            IF(GamRL.GE.0.7) THEN
            EPSpl=EPSpl0
            Suln=Suln0
            Ncyc=Ncyc0
            END IF
            END IF
          GO TO 10
          ELSE IF(Jcon0.EQ.4) THEN
          Jcon=4
          Jcon0=3
            IF(Sigrl.EQ.0.0) THEN
            IF(GamRL.GE.0.7) THEN
            EPSpl=EPSpl0
            Suln=Suln0
            Ncyc=Ncyc0
            END IF
            END IF
          GO TO 5
          ELSE
          Jcon=1
          Jcon0=3
          GO TO 20
          END IF
C
C  Otherwise, still reloading along the previous unloading path
        ELSE
        Sig1=Sigul*((EPS1-EPSpl)/(EPSul-EPSpl))**2.0
        YMtan=2*Sigul*(EPS1-EPSpl)/((EPSul-EPSpl)**2.0)
        END IF
C
C *** RELOADING (NORMALIZED STRAIN IS LARGER THAN 1) (Jcon=7) *****
C *** Hysteresis has exceeded the unloaded strain from the envelope curve
C  Jcon0 in not renewed here. (Jcon0=3)
      ELSE IF(Jcon.EQ.7) THEN
      GO TO 10
C
C *** VIRGIN CONCRETE (Not Loaded in Compression) (Jcon=8) ***********
      ELSE
        IF(EPS1.GE.0.0) THEN
        Sig1=0.0
        YMtan=0.0
        ELSE
        Jcon=1
        Jcon0=8
        GO TO 20
        END IF
C --------------------------------------------------------------------
      END IF
      GO TO 30
C
C ********************************************************************
C *** Reloading Path (Linear Function, 0.2< Normalized Strain <1)
    5 CONTINUE
      IF(EPS1.LT.EPSule) THEN
      Jcon=7
      GO TO 10
      ELSE
      Sig1=YMrl*(EPS1-EPSule)+Suln
      YMtan=YMrl
      END IF
      GO TO 30
C
C ********************************************************************
C *** When Jcon=7
   10 CONTINUE
C
C *** Stress on the reloading path
      Sign0=YMrl*(EPS1-EPSule)+Suln
      YMtn0=YMrl
C
C *** Stress on the envelope curve (Mander Model)
      EPSxx=EPS1/EPScc
      AAA=EPSxx**r
      Senv=Sigcc*EPSxx*r/(r-1.0+AAA)
      YMtenv=Sigcc/EPScc*r*(r-1.0)*(1.0-AAA)/(r-1.0+AAA)**2.0
C
C *** If hysteresis has reached the envelope curve or not
      IF(Sign0.GT.Senv) THEN
      Sig1=Sign0
      YMtan=YMtn0
      ELSE
      Jcon=1
      Jcon0=7
      Sig1=Senv
      YMtan=YMtenv
      END IF
      GO TO 30
C
C *** ENVELOPE CURVE (Mander Model)
   20 CONTINUE
      EPSxx=EPS1/EPScc
      AAA=EPSxx**r
      Sig1=Sigcc*EPSxx*r/(r-1.0+AAA)
      YMtan=Sigcc/EPScc*r*(r-1.0)*(1.0-AAA)/(r-1.0+AAA)**2.0
      GO TO 30
C --------------------------------------------------------------------
C
   30 CONTINUE
      RETURN
      END
C ====================================================================
C
C ====================================================================
C ----------------------- SUBROUTINE PROGRAM -------------------------
C ======================= UNLOADING (TENSION) ========================
C ====================================================================
C
      SUBROUTINE TENSI14 (EPS1,Sig1,YMtan,Jcon,EPS0,Sig0,
     1                    EPSule,Sigule,EPSpl,Suln,EPSul,Sigul,
     1                    EPSrl,Sigrl,EPSpl0,Suln0,GamRL,
     1                    Ncyc,Jcon0,Ncyc0)
      IMPLICIT REAL*8(A-H,O-Z)
C
C *** UNLOADING FROM THE ENVELOPE CURVE (Jcon=1) *********************
      IF(Jcon.EQ.1) THEN
      Ncyc=1
      Jcon=3
      Jcon0=1
      GO TO 10
C
C *** ZERO STRESS (Jcon=2) *******************************************
      ELSE IF(Jcon.EQ.2) THEN
      Sig1=0.0
      YMtan=0.0
C
C *** ON UNLOADING PATH (Jcon=3) *************************************
      ELSE IF(Jcon.EQ.3) THEN
        IF(EPS1.GE.EPSpl) THEN
        Jcon=2
        Jcon0=3
        Sig1=0.0
        YMtan=0.0
        ELSE
        Sig1=Sigul*((EPS1-EPSpl)/(EPSul-EPSpl))**2.0
        YMtan=2.0*Sigul*(EPS1-EPSpl)/((EPSul-EPSpl)**2.0)
        END IF
C
C *** UNLOADING FROM RELOADING PATH (Jcon=4) *************************
      ELSE  IF(Jcon.EQ.4) THEN
      Jcon=3
      Jcon0=4
      GO TO 20
C
C *** ON UNLOADING PATH, 0< Normalized Strain <0.2 (Jcon=5) **********
      ELSE IF(Jcon.EQ.5) THEN
        IF(EPS1.GE.EPSpl) THEN
        Jcon=2
        Sig1=0.0
        YMtan=0.0
        ELSE
        Sig1=2.5*Suln0*((EPS1-EPSpl)/(EPSule-EPSpl))**2.0
        YMtan=5.0*Suln0*(EPS1-EPSpl)/((EPSule-EPSpl)**2.0)
        END IF
C
C *** UNLOADING FROM RELOADING ALONG THE PREVIOUS UNLOADING PATH (Jcon=6)
      ELSE IF(Jcon.EQ.6) THEN
C
        IF(EPS1.GE.EPSpl) THEN
        Jcon0=3
        Jcon=2
        Sig1=0.0
        YMtan=0.0
        ELSE
C  Jcon0 in not renewed here.
        Jcon=3
        Sig1=Sigul*((EPS1-EPSpl)/(EPSul-EPSpl))**2.0
        YMtan=2.0*Sigul*(EPS1-EPSpl)/((EPSul-EPSpl)**2.0)
        END IF
C
C *** UNLOADING FROM WHERE NORMALIZED STRAIN IS LARGER THAN 1 (Jcon=7)
      ELSE IF(Jcon.EQ.7) THEN
C *** (EPSule, Sigule) is renewed.
      EPSule=EPS0
      Sigule=Sig0
C
      Suln=Sig0
      Jcon=3
      Jcon0=7
      GO TO 20
C
C *** VIRGIN CONCRETE (Not Loaded in Compression) (Jcon=8) ***********
      ELSE
      Sig1=0.0
      YMtan=0.0
C --------------------------------------------------------------------
      END IF
      GO TO 30
C
C ********************************************************************
C *** UNLOADING FROM THE ENVELOPE CURVE
   10 CONTINUE
C *** Initialize the data about unloading point
      EPSule=EPS0
      Sigule=Sig0
C
      EPSul=EPSule
      Sigul=Sigule
C
      EPSrl=EPSule
      Sigrl=Sigule
C
C *** Prediction of the Plastic Strain from Eq.(9)
        IF(EPSule.GT.-1000.0e-6) THEN
        EPSpl=0.0
        ELSE IF(EPSule.GT.-3500.0e-6) THEN
        EPSpl=0.43*(EPSule+0.001)
        ELSE
        EPSpl=0.94*(EPSule+0.00235)
        END IF
C
        EPSpl0=EPSpl
C
C
C *** Prediction of the Unloaded Stress from Eqs.(18) and(20)
        IF(EPSule.GT.-1000.0e-6) THEN
        Suln=Sigule*1.0
        ELSE IF(EPSule.GT.-3500.0e-6) THEN
        Suln=Sigule*(1.0+32.0*(EPSule+0.001))
        ELSE
        Suln=Sigule*0.92
        END IF
C
        Suln0=Sigul
        Ncyc0=Ncyc
        GamRL=0.0
C
C *** Already reached the plastic strain
        IF(EPS1.GE.EPSpl) THEN
        Jcon=2
        Jcon0=3
        Sig1=0.0
        YMtan=0.0
C *** Otherwise, on unloading path
        ELSE
        Sig1=Sigul*((EPS1-EPSpl)/(EPSul-EPSpl))**2.0
        YMtan=2.0*Sigul*(EPS1-EPSpl)/((EPSul-EPSpl)**2.0)
        END IF
      GO TO 30
C
C ********************************************************************
C *** UNLOADING FROM RELOADING PATH
   20 CONTINUE
C *** Initialize the data about unloaded point
      EPSul=EPS0
      Sigul=Sig0
      GamRL=(EPS0-EPSpl)/(EPSule-EPSpl)
C
C ======= If it satisfies Eq. (6), the number of cycles is counted.
        IF(Sigrl.EQ.0.0) THEN
        IF(GamRL.GE.0.7) THEN
        Ncyc0=Ncyc
        Ncyc=Ncyc+1
C
          IF(GamRL.GT.1.0) THEN
          GamRL=1.0
          END IF
C
C *** Prediction of the Plastic Strain from Eqs.(11), (12) and (13)
        EPSpl0=EPSpl
          IF(Ncyc.LE.2) THEN
          Gamma=0.945+0.2*(1.0-GamRL)
          ELSE
          Gamma=0.965+0.005*(Ncyc-3)+0.2*(1.0-GamRL)
          END IF
C
          IF(Gamma.GT.1.0) THEN
          Gamma=1.0
          END IF
C
          IF(EPSule.GT.-1000.0e-6) THEN
          Gamma=1.0
          END IF
C
        EPSpl=(1.0-Gamma)*EPSule+Gamma*EPSpl
C
C
C *** Prediction of the Unloaded Stress from Eqs.(18) ~ (21)
          Suln0=Suln
C
          IF(Ncyc.LE.2) THEN
            IF(EPSule.GT.-3500.0e-6) THEN
            aBeta=1.0+(42.0-10.0*Ncyc)*(EPSule+0.001)
            ELSE
            aBeta=0.945
            END IF
          ELSE
            IF(EPSule.GT.-3500.0e-6) THEN
            aBeta=1.0+(20.0-2.0*Ncyc)*(EPSule+0.001)
            ELSE
            aBeta=0.965+0.005*(Ncyc-3)
            END IF
          END IF
          Beta=aBeta+0.2*(1.0-GamRL)
C
          IF(Beta.GT.1.0) THEN
          Beta=1.0
          END IF
C
          Suln=Suln*Beta
C
C ======= If it does not satisfy Eq. (6),
C         and if the partial reloading ratio is less than 0.7, then
        ELSE
        Suln=Suln
        EPSpl=EPSpl
        END IF
        GO TO 22
C
C         If the reloading is not from zero stress, then
        ELSE
        Suln=Suln
        EPSpl=EPSpl
        END IF
C
   22   IF(EPS1.GE.EPSpl) THEN
        Jcon=2
        Jcon0=3
        Sig1=0.0
        YMtan=0.0
        ELSE
        Sig1=Sigul*((EPS1-EPSpl)/(EPSul-EPSpl))**2.0
        YMtan=2.0*Sigul*(EPS1-EPSpl)/((EPSul-EPSpl)**2.0)
        END IF
      GO TO 30
C --------------------------------------------------------------------
   30 CONTINUE
      RETURN
      END
C ====================================================================
C ====================================================================
C ====================================================================

max_stars_repo_name: kuanshi/ductile-fracture

max_issues_repo_name: kuanshi/ductile-fracture

max_forks_repo_name: kuanshi/ductile-fracture


text:
!
!  Copyright 2017-2019 SALMON developers
!
!  Licensed under the Apache License, Version 2.0 (the "License");
!  you may not use this file except in compliance with the License.
!  You may obtain a copy of the License at
!
!      http://www.apache.org/licenses/LICENSE-2.0
!
!  Unless required by applicable law or agreed to in writing, software
!  distributed under the License is distributed on an "AS IS" BASIS,
!  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
!  See the License for the specific language governing permissions and
!  limitations under the License.
!
!=======================================================================
!=======================================================================

subroutine set_vonf_sd
!$ use omp_lib
  use salmon_parallel, only: nproc_id_global
  use salmon_communication, only: comm_is_root, comm_summation
  use scf_data
  implicit none
  integer :: i
  integer :: ix,iy,iz,iix,iiy,iiz
  integer :: max_icell
  real(8) :: rr

  if(iperiodic==0)then
    max_icell=0
  else if(iperiodic==3)then
    max_icell=2
  end if

  vonf_sd=0.d0
  eonf_sd=0.d0

  do i=1,nump
    do iiz=-max_icell,max_icell
    do iiy=-max_icell,max_icell
    do iix=-max_icell,max_icell
      do iz=mg_sta(3),mg_end(3)
      do iy=mg_sta(2),mg_end(2)
      do ix=mg_sta(1),mg_end(1)
        rr=sqrt((gridcoo(ix,1)-(coop(1,i)+dble(iix*lg_num(1))*Hgs(1)))**2 &
               +(gridcoo(iy,2)-(coop(2,i)+dble(iiy*lg_num(2))*Hgs(2)))**2 &
               +(gridcoo(iz,3)-(coop(3,i)+dble(iiz*lg_num(3))*Hgs(3)))**2)
        if(rr>radp_diele)then
          vonf_sd(ix,iy,iz)=vonf_sd(ix,iy,iz)  &
                        -(vecp(1,i)*(gridcoo(ix,1)-(coop(1,i)+dble(iix*lg_num(1))*Hgs(1))) &
                         +vecp(2,i)*(gridcoo(iy,2)-(coop(2,i)+dble(iiy*lg_num(2))*Hgs(2))) &
                         +vecp(3,i)*(gridcoo(iz,3)-(coop(3,i)+dble(iiz*lg_num(3))*Hgs(3))))/rr**3
          eonf_sd(1,ix,iy,iz)=eonf_sd(1,ix,iy,iz)  &
                          +(3.d0*(vecp(1,i)*(gridcoo(ix,1)-(coop(1,i)+dble(iix*lg_num(1))*Hgs(1)))/rr &
                                 +vecp(2,i)*(gridcoo(iy,2)-(coop(2,i)+dble(iiy*lg_num(2))*Hgs(2)))/rr &
                                 +vecp(3,i)*(gridcoo(iz,3)-(coop(3,i)+dble(iiz*lg_num(3))*Hgs(3)))/rr)* &
                              (gridcoo(ix,1)-(coop(1,i)+dble(iix*lg_num(1))*Hgs(1)))/rr &
                        -vecp(1,i)) /rr**3
          eonf_sd(2,ix,iy,iz)=eonf_sd(2,ix,iy,iz)  &
                          +(3.d0*(vecp(1,i)*(gridcoo(ix,1)-(coop(1,i)+dble(iix*lg_num(1))*Hgs(1)))/rr &
                                 +vecp(2,i)*(gridcoo(iy,2)-(coop(2,i)+dble(iiy*lg_num(2))*Hgs(2)))/rr &
                                 +vecp(3,i)*(gridcoo(iz,3)-(coop(3,i)+dble(iiz*lg_num(3))*Hgs(3)))/rr)* &
                              (gridcoo(iy,2)-(coop(2,i)+dble(iiy*lg_num(2))*Hgs(2)))/rr &
                        -vecp(2,i)) /rr**3
          eonf_sd(3,ix,iy,iz)=eonf_sd(3,ix,iy,iz)  &
                          +(3.d0*(vecp(1,i)*(gridcoo(ix,1)-(coop(1,i)+dble(iix*lg_num(1))*Hgs(1)))/rr &
                                 +vecp(2,i)*(gridcoo(iy,2)-(coop(2,i)+dble(iiy*lg_num(2))*Hgs(2)))/rr &
                                 +vecp(3,i)*(gridcoo(iz,3)-(coop(3,i)+dble(iiz*lg_num(3))*Hgs(3)))/rr)* &
                              (gridcoo(iz,3)-(coop(3,i)+dble(iiz*lg_num(3))*Hgs(3)))/rr &
                          -vecp(3,i)) /rr**3
        else
          vonf_sd(ix,iy,iz)=vonf_sd(ix,iy,iz)  &
                        -(vecp(1,i)*(gridcoo(ix,1)-(coop(1,i)+dble(iix*lg_num(1))*Hgs(1))) &
                         +vecp(2,i)*(gridcoo(iy,2)-(coop(2,i)+dble(iiy*lg_num(2))*Hgs(2))) &
                         +vecp(3,i)*(gridcoo(iz,3)-(coop(3,i)+dble(iiz*lg_num(3))*Hgs(3))))/radp_diele**3
          eonf_sd(1,ix,iy,iz)=eonf_sd(1,ix,iy,iz)  &
                         +vecp(1,i)/radp_diele**3
          eonf_sd(2,ix,iy,iz)=eonf_sd(2,ix,iy,iz)  &
                         +vecp(2,i)/radp_diele**3
          eonf_sd(3,ix,iy,iz)=eonf_sd(3,ix,iy,iz)  &
                         +vecp(3,i)/radp_diele**3
        end if
      end do
      end do
      end do
    end do
    end do
    end do
  end do

  return

end subroutine set_vonf_sd

max_stars_repo_name: shunsuke-sato/SALMON

max_issues_repo_name: shunsuke-sato/SALMON

max_forks_repo_name: shunsuke-sato/SALMON


text:
      SUBROUTINE STRUCTM(Xi,Scale2,Uv,Dv,Us,Ds,Ss,Cs,Bs,Ts,Gl)
      IMPLICIT NONE
      DOUBLE PRECISION Bs , Cs , Ds , Dv , Gl , Scale2 , Ss , Ts , Us ,
     &                 Uv , Xi
      END SUBROUTINE

ext: f

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET