diff --git a/lib/neuronet/ryo_force_NN.F90 b/lib/neuronet/ryo_force_NN.F90
index b82c85f98..ecd83937f 100644
--- a/lib/neuronet/ryo_force_NN.F90
+++ b/lib/neuronet/ryo_force_NN.F90
@@ -1,530 +1,530 @@
 module NN
      !-----------------------------------------------------------------------
      !                        Time-stamp: <2016-09-06 16:34:13 Ryo KOBAYASHI>
      !-----------------------------------------------------------------------
      !  Parallel implementation of neural-network potential with 1 hidden
      !  layer. It is available for plural number of species.
      !-----------------------------------------------------------------------
      !.....cutoff region width ratio to rc
      real(8):: rcw = 0.9d0
 
 
 contains
      subroutine force_NN(namax,natm,tag,ra,nnmax,force,vatom,h,hi,tcom &
           ,nb,nbmax,lsb,nex,lsrc,myparity,nn,sv,rc,lspr &
           ,mpi_world,myid,epi,epot,vflag_atom, itype, cnst, iaddr2, iaddr3, &
           nsp, nhl, max_ncnst, nl, nlmax, hl1, hl2, gsf, dgsf, nal, nnl, wgt11, &
           wgt12, wgt21, wgt22, wgt23, rc3)
           ! namax = nmax
           ! natm = inum
           ! tag = type
           ! ra = x ! no need to tranpose
           ! nnmax = nnmax
           ! force = fort_f ! no need to tranpose
           ! vatom = fort_vatom
           ! h = h
           ! hi = h
           ! tcom = dummy
           ! nb = nghost
           ! nbmax = idummy
           ! lsb = idummy
           ! nex = idummy
           ! lsrc = idummy
           ! myparity = idummy
           ! nn = idummy
           ! sv = dummy
           ! rc = rcin
           ! lspr = fort_neighbor
           ! mpi_world = idummy
           ! myid = idummy
           ! epi = fort_eatom
           ! epot = eng_ndwl
           ! vflag_atom = vflag_atom
 
           implicit none
           !include "mpif.h"
           !include "./params_unit.h"
           integer,intent(in):: nhl(0:nlmax+1), itype(nhl(0)), nl, nlmax, max_ncnst, nsp, nal, nnl
           real(8),intent(in):: cnst(max_ncnst,nhl(0))
           integer,intent(in):: iaddr2(2,nsp,nsp),iaddr3(2,nsp,nsp,nsp)
           real(8),intent(in):: wgt11(nhl(0), nhl(1)),wgt12(nhl(1))
           real(8),intent(in):: wgt21(nhl(0),nhl(1)),wgt22(nhl(1),nhl(2)), &
                wgt23(nhl(2))
 
           integer,intent(in):: namax,natm,nnmax
           integer,intent(in):: nb,nbmax,lsb(0:nbmax,6),lsrc(6),myparity(3) &
                ,nn(6),mpi_world,myid,lspr(0:nnmax,namax),nex(3),tag(namax)
           real(8),intent(in):: ra(3,namax),h(3,3),hi(3,3),sv(3,6)
           real(8),intent(inout):: tcom,rc
           real(8),intent(out):: force(3,namax),epi(namax),epot,vatom(6,namax)
           logical:: vflag_atom
 
           !.....local
           integer:: i,j,k,l,m,n,ixyz,jxyz,is,js,ks,ierr,nbl,ia,ja,nexp,isf &
                ,icoeff,ihl0,ihl1,ihl2,jj,jsf
           real(8):: rcin,b_na,at(3),epotl,wgt,hl1i,hl2i,tmp2,tmp1,tmp,tmp3(3)
           real(8),intent(in):: rc3
           real(8),intent(inout) :: hl1(nhl(1),nal),hl2(nhl(2),nal)
           real(8),intent(inout) :: gsf(nhl(0),nal),dgsf(3,nhl(0),0:nnl,nal)
           !    real(8),allocatable:: aml(:,:,:,:),bml(:,:,:,:)
           !.....1st call
           logical,save:: l1st=.true.
 
 
           !.....first, calculate all the symmetry functions
           call eval_sf(nhl(0),namax,natm,nb,nnmax,h,tag,ra &
                ,lspr,rc,rc3, itype, cnst, nhl, nlmax, max_ncnst &
                ,iaddr2, iaddr3, nsp, gsf, dgsf, nal, nnl)
 
           !.....2nd, calculate the node values by summing contributions from
           !.....  symmetry functions
           if( nl.eq.1 ) then
              hl1(1:nhl(1),1:natm)= 0d0
              do ia=1,natm
                 !.....debug
                 do ihl1=1,nhl(1)
                    tmp= 0d0
                    do ihl0=1,nhl(0)
                       tmp= tmp +wgt11(ihl0,ihl1) *gsf(ihl0,ia)
                    enddo
                    hl1(ihl1,ia)= sigmoid(tmp)
                 enddo
              enddo
           else if( nl.eq.2 ) then
              hl1(1:nhl(1),1:natm)= 0d0
              hl2(1:nhl(2),1:natm)= 0d0
              do ia=1,natm
                 do ihl1=1,nhl(1)
                    tmp= 0d0
                    do ihl0=1,nhl(0)
                       tmp= tmp +wgt21(ihl0,ihl1) *gsf(ihl0,ia)
                    enddo
                    hl1(ihl1,ia)= sigmoid(tmp)
                 enddo
                 do ihl2=1,nhl(2)
                    tmp= 0d0
                    do ihl1=1,nhl(1)
                       tmp= tmp +wgt22(ihl1,ihl2) *(hl1(ihl1,ia)-0.5d0)
                    enddo
                    hl2(ihl2,ia)= sigmoid(tmp)
                 enddo
              enddo
 
           endif
 
           !.....then calculate the energy of atom by summing up the node values
           epotl= 0d0
           if( nl.eq.1 ) then
              do ia=1,natm
                 epi(ia)= 0d0
                 do ihl1=1,nhl(1)
                    epi(ia)= epi(ia) +wgt12(ihl1) *(hl1(ihl1,ia)-0.5d0)
                 enddo
                 epotl=epotl +epi(ia)
              enddo
           else if( nl.eq.2 ) then
              do ia=1,natm
                 epi(ia)= 0d0
                 do ihl2=1,nhl(2)
                    epi(ia)= epi(ia) +wgt23(ihl2) *(hl2(ihl2,ia)-0.5d0)
                 enddo
                 epotl=epotl +epi(ia)
              enddo
           endif
 
           !.....sum up for forces
           force(1:3,1:natm+nb)= 0d0
           if( nl.eq.1 ) then
              do ia=1,natm
                 do ihl1=1,nhl(1)
                    hl1i= hl1(ihl1,ia)
                    tmp= wgt12(ihl1)*hl1i*(1d0-hl1i)
                    do jj=1,lspr(0,ia)
                       ja= lspr(jj,ia)
                       do ihl0=1,nhl(0)
                          force(1:3,ja)=force(1:3,ja) &
                               -tmp*wgt11(ihl0,ihl1)*dgsf(1:3,ihl0,jj,ia)
                       enddo
                    enddo
                    !.....atom ia
                    do ihl0= 1,nhl(0)
                       force(1:3,ia)=force(1:3,ia) &
                            -tmp*wgt11(ihl0,ihl1)*dgsf(1:3,ihl0,0,ia)
                    enddo
                 enddo
              enddo
           else if( nl.eq.2 ) then
              do ia=1,natm
                 do ihl2=1,nhl(2)
                    hl2i= hl2(ihl2,ia)
                    tmp2= wgt23(ihl2) *hl2i*(1d0-hl2i)
                    do ihl1=1,nhl(1)
                       hl1i= hl1(ihl1,ia)
                       tmp1= wgt22(ihl1,ihl2) *hl1i*(1d0-hl1i)
                       do jj=1,lspr(0,ia)
                          ja= lspr(jj,ia)
                          do ihl0=1,nhl(0)
                             force(1:3,ja)=force(1:3,ja) &
                                  -tmp2 *tmp1 &
                                  *wgt21(ihl0,ihl1)*dgsf(1:3,ihl0,jj,ia)
                          enddo
                       enddo
                       !.....atom ia
                       do ihl0= 1,nhl(0)
                          force(1:3,ia)=force(1:3,ia) &
                               -tmp2*tmp1*wgt21(ihl0,ihl1)*dgsf(1:3,ihl0,0,ia)
                       enddo
                    enddo
                 enddo
              enddo
           endif
 
           !          call copy_dba_bk(tcom,namax,natm,nbmax,nb,lsb,nex,lsrc,myparity &
           !               ,nn,mpi_world,force,3)
 !!$    if( myid.ge.0 ) then
 !!$      call copy_dba_bk(tcom,namax,natm,nbmax,nb,lsb,lsrc,myparity &
 !!$           ,nn,mpi_world,force,3)
 !!$    else
 !!$      call reduce_dba_bk(natm,namax,tag,force,3)
 !!$    endif
 
 
           if( vflag_atom ) then
              call compute_stress(namax,natm,tag,ra,nnmax,vatom,h &
                   ,tcom,nb,nbmax,lsb,nex,lsrc,myparity,nn,rc,lspr &
                   ,mpi_world,myid, nl, wgt11, wgt12, wgt21, wgt22 &
                   , wgt23, nlmax, nhl, hl1, hl2, nal, dgsf, nnl)
           endif
 
           epot= epotl
 
           return
      end subroutine force_NN
      !=======================================================================
      subroutine eval_sf(nsf,namax,natm,nb,nnmax,h,tag,ra,lspr,rc,rc3, itype,&
           cnst, nhl, nlmax, max_ncnst, iaddr2, iaddr3, nsp, gsf, dgsf, nal, nnl)
           !
           !  Evaluate symmetry functions and derivatives for multi-species system.
           !
           implicit none
           integer,intent(in):: nsf,namax,natm,nb,nnmax,lspr(0:nnmax,namax),&
                tag(namax), nlmax, nhl(0:nlmax+1), itype(nhl(0)), max_ncnst,&
                nsp, nal, nnl
           integer,intent(in):: iaddr2(2,nsp,nsp),iaddr3(2,nsp,nsp,nsp)
           real(8),intent(in):: h(3,3),ra(3,namax),rc,rc3
           !    real(8),intent(out):: gsf(nsf,nal),dgsf(3,nsf,0:nnl,nal)
           real(8),intent(in):: cnst(max_ncnst,nhl(0))
           integer:: isf,isfc,ia,jj,ja,kk,ka,is,js,ks,isfc1,isfc2
           real(8):: xi(3),xj(3),xij(3),rij(3),dij,fcij,eta,rs,texp,driji(3), &
                dfcij,drijj(3),dgdr,xk(3),xik(3),rik(3),dik,fcik,dfcik, &
                driki(3),drikk(3),almbd,spijk,cs,t1,t2,dgdij,dgdik,dgcs, &
                dcsdj(3),dcsdk(3),dcsdi(3),tcos,tpoly,a1,a2,tmorse
 
           real(8),external:: sprod
           real(8),intent(out) :: gsf(nhl(0),nal), dgsf(3,nhl(0),0:nnl,nal)
 
           gsf(1:nsf,1:nal)= 0d0
           dgsf(1:3,1:nsf,0:nnl,1:nal)= 0d0
           do ia=1,natm
-             xi(1:3)= ra(1:3,ia)
-             is= tag(ia)
+             xi(1:3)= ra(1:3,ia) ! pos atom_i
+             is= tag(ia)         ! type atom_i
              do jj=1,lspr(0,ia)
-                ja= lspr(jj,ia)
+                ja= lspr(jj,ia)  ! neighbor index j
                 if( ja.eq.ia ) cycle
-                xj(1:3)= ra(1:3,ja)
-                xij(1:3)= xj(1:3)-xi(1:3)
-                rij(1:3)= h(1:3,1)*xij(1) +h(1:3,2)*xij(2) +h(1:3,3)*xij(3)
-                dij= sqrt(rij(1)**2 +rij(2)**2 +rij(3)**2)
+                xj(1:3)= ra(1:3,ja)  ! pos atom_j
+                xij(1:3)= xj(1:3)-xi(1:3) ! vec from atom i to j
+                rij(1:3)= h(1:3,1)*xij(1) +h(1:3,2)*xij(2) +h(1:3,3)*xij(3) ! unscale 
+                dij= sqrt(rij(1)**2 +rij(2)**2 +rij(3)**2) ! distance ij
                 if( dij.ge.rc ) cycle
-                js= tag(ja)
+                js= tag(ja) ! type atom j
                 isfc=0
-                driji(1:3)= -rij(1:3)/dij
-                drijj(1:3)= -driji(1:3)
-                fcij= fc(dij,rc)
-                dfcij= dfc(dij,rc)
+                driji(1:3)= -rij(1:3)/dij !unit vec i->j
+                drijj(1:3)= -driji(1:3)   !unit vec j->i
+                fcij= fc(dij,rc)          !eq (4)
+                dfcij= dfc(dij,rc)        !eq (5)
                 do isf=iaddr2(1,is,js),iaddr2(2,is,js)
 !!$          isfc= isfc+1
 !!$          isf= (icmb2(is,js)-1)*nsfc1 +isfc1
                    if( itype(isf).eq.1 ) then ! Gaussian
                       eta= cnst(1,isf)
                       rs=  cnst(2,isf)
                       !.....function value
                       texp= exp(-eta*(dij-rs)**2)
                       gsf(isf,ia)= gsf(isf,ia) +texp*fcij
                       !.....derivative
                       ! dgsf(ixyz,isf,jj,ia): derivative of isf-th basis of atom-ia
                       ! by ixyz coordinate of atom-jj.
                       ! jj=0 means derivative by atom-ia.
                       dgdr= -2d0*eta*(dij-rs)*texp*fcij +texp*dfcij
                       dgsf(1:3,isf,0,ia)= dgsf(1:3,isf,0,ia) +driji(1:3)*dgdr
                       dgsf(1:3,isf,jj,ia)= dgsf(1:3,isf,jj,ia) +drijj(1:3)*dgdr
                    else if( itype(isf).eq.2 ) then ! cosine
                       a1= cnst(1,isf)
                       !.....func value
                       tcos= (1d0+cos(dij*a1))
                       gsf(isf,ia)= gsf(isf,ia) +tcos*fcij
                       !.....derivative
                       dgdr= -a1*sin(dij*a1)*fcij +tcos*dfc(dij,rc)
                       dgsf(1:3,isf,0,ia)= dgsf(1:3,isf,0,ia) +driji(1:3)*dgdr
                       dgsf(1:3,isf,jj,ia)= dgsf(1:3,isf,jj,ia) +drijj(1:3)*dgdr
                    else if( itype(isf).eq.3 ) then ! polynomial
                       a1= cnst(1,isf)
                       !.....func value
                       tpoly= 1d0*dij**(-a1)
                       gsf(isf,ia)= gsf(isf,ia) +tpoly*fcij
                       !.....derivative
                       dgdr= -a1*dij**(-a1-1d0)*fcij +tpoly*dfc(dij,rc)
                       dgsf(1:3,isf,0,ia)= dgsf(1:3,isf,0,ia) +driji(1:3)*dgdr
                       dgsf(1:3,isf,jj,ia)= dgsf(1:3,isf,jj,ia) +drijj(1:3)*dgdr
                    else if( itype(isf).eq.4 ) then ! Morse-type
                       a1= cnst(1,isf)
                       a2= cnst(2,isf)
                       !.....func value
                       texp= exp(-a1*(dij-a2))
                       tmorse= ((1d0-texp)**2 -1d0)
                       gsf(isf,ia)= gsf(isf,ia) +tmorse*fcij
                       !.....derivative
                       dgdr= 2d0*a1*(1d0-texp)*texp*fcij +tmorse*dfcij
                       dgsf(1:3,isf,0,ia)= dgsf(1:3,isf,0,ia) +driji(1:3)*dgdr
                       dgsf(1:3,isf,jj,ia)= dgsf(1:3,isf,jj,ia) +drijj(1:3)*dgdr
                    endif
                 enddo
 
 !!$        fcij= fc(dij,rc)
 !!$        dfcij= dfc(dij,rc)
 !!$        driji(1:3)= -rij(1:3)/dij
 !!$        drijj(1:3)= -driji(1:3)
                 if( dij.gt.rc3 ) cycle
                 do kk=1,lspr(0,ia)
                    ka= lspr(kk,ia)
                    ks= tag(ka)
                    if( iaddr3(1,is,js,ks).lt.0 ) cycle
                    if( ka.eq.ia .or. ka.le.ja ) cycle
                    xk(1:3)= ra(1:3,ka)
                    xik(1:3)= xk(1:3)-xi(1:3)
                    rik(1:3)= h(1:3,1)*xik(1) +h(1:3,2)*xik(2) +h(1:3,3)*xik(3)
                    dik= sqrt(rik(1)**2 +rik(2)**2 +rik(3)**2)
                    if( dik.ge.rc3 ) cycle
                    do isf=iaddr3(1,is,js,ks),iaddr3(2,is,js,ks)
 !!$            isf= nsfc1*nc1 +(icmb3(is,js,ks)-1)*nsfc2 +isfc2
                       almbd= cnst(1,isf)
                       t2= (abs(almbd)+1d0)**2
                       fcik= fc(dik,rc)
                       dfcik= dfc(dik,rc)
                       driki(1:3)= -rik(1:3)/dik
                       drikk(1:3)= -driki(1:3)
                       !.....function value
                       spijk= rij(1)*rik(1) +rij(2)*rik(2) +rij(3)*rik(3)
                       cs= spijk/dij/dik
                       t1= (almbd +cs)**2
                       gsf(isf,ia)= gsf(isf,ia) +t1/t2 *fcij*fcik
                       !.....derivative
                       dgdij= dfcij *fcik *t1/t2
                       dgdik= fcij *dfcik *t1/t2
                       dgsf(1:3,isf,0,ia)= dgsf(1:3,isf,0,ia) &
                            +dgdij*driji(1:3) +dgdik*driki(1:3)
                       dgsf(1:3,isf,jj,ia)= dgsf(1:3,isf,jj,ia) +dgdij*drijj(1:3)
                       dgsf(1:3,isf,kk,ia)= dgsf(1:3,isf,kk,ia) +dgdik*drikk(1:3)
                       dgcs= 2d0*(almbd+cs)/t2 *fcij*fcik
                       dcsdj(1:3)= rik(1:3)/dij/dik -rij(1:3)*spijk/dij**3/dik
                       dcsdk(1:3)= rij(1:3)/dij/dik -rik(1:3)*spijk/dik**3/dij
                       dcsdi(1:3)= -dcsdj(1:3) -dcsdk(1:3)
                       dgsf(1:3,isf,0,ia)= dgsf(1:3,isf,0,ia) +dgcs*dcsdi(1:3)
                       dgsf(1:3,isf,jj,ia)= dgsf(1:3,isf,jj,ia) +dgcs*dcsdj(1:3)
                       dgsf(1:3,isf,kk,ia)= dgsf(1:3,isf,kk,ia) +dgcs*dcsdk(1:3)
                    enddo
                 enddo
              enddo
           enddo
 
      end subroutine eval_sf
      !=======================================================================
      function fc(r,rc)
           implicit none
           real(8),intent(in):: r,rc
           real(8):: fc,rs
           real(8),parameter:: pi= 3.14159265358979d0
           rs= rc*rcw
           if( r.le.rs ) then
              fc= 1d0
           else if( r.gt.rs .and. r.le.rc ) then
              fc= 0.5d0 *(cos((r-rs)/(rc-rs)*pi)+1d0)
           else
              fc= 0d0
           endif
           return
      end function fc
      !=======================================================================
      function dfc(r,rc)
           implicit none
           real(8),intent(in):: r,rc
           real(8):: dfc,rs
           real(8),parameter:: pi= 3.14159265358979d0
           rs= rc*rcw
           if( r.le.rs ) then
              dfc= 0d0
           else if( r.gt.rs .and. r.le.rc ) then
              dfc= -pi/2/(rc-rs) *sin((r-rs)/(rc-rs)*pi)
           else
              dfc= 0d0
           endif
           return
      end function dfc
      !=======================================================================
      function sigmoid(x)
           implicit none
           real(8),intent(in):: x
           real(8):: sigmoid
 
           sigmoid= 1d0/(1d0 +exp(-x))
           return
      end function sigmoid
      !=======================================================================
      function dsigmoid(x)
           implicit none
           real(8),intent(in):: x
           real(8):: dsigmoid,sx
 
           sx= sigmoid(x)
           !    dsigmoid= -exp(-x)/(1d0+exp(-x))**2
           dsigmoid= sx*(1d0-sx)
           return
      end function dsigmoid
      !=======================================================================
      function factorial(n,m)
           !  compute factorial of n, m-times.
           implicit none
           integer,intent(in):: n,m
           real(8):: factorial
 
           integer:: i
 
           factorial= 1
           do i=0,m-1
              factorial= factorial*(n-i)
           enddo
           return
      end function factorial
      !=======================================================================
      subroutine compute_stress(namax,natm,tag,ra,nnmax,vatom,h &
           ,tcom,nb,nbmax,lsb,nex,lsrc,myparity,nn,rc,lspr &
           ,mpi_world,myid, nl, wgt11, wgt12, wgt21, wgt22, wgt23, nlmax, nhl,&
           hl1, hl2, nal, dgsf, nnl)
           implicit none
           integer,intent(in):: namax,natm,nnmax,nb,nbmax,lsb(0:nbmax,6)&
                ,lsrc(6),myparity(3),nn(6),mpi_world,myid,lspr(0:nnmax,namax)&
                ,nex(3),tag(namax), nl, nlmax, nhl(0:nlmax+1), nal, nnl
           real(8),intent(in):: ra(3,namax),h(3,3),rc, dgsf(3,nhl(0),0:nnl,nal)
           real(8),intent(inout):: tcom
           real(8),intent(out):: vatom(6,namax)
           real(8),intent(in):: wgt11(nhl(0), nhl(1)),wgt12(nhl(1))
           real(8),intent(in):: wgt21(nhl(0),nhl(1)),wgt22(nhl(1),nhl(2)), &
                wgt23(nhl(2))
           real(8),intent(in) :: hl1(nhl(1),nal)
           real(8),intent(in):: hl2(nhl(2),nal)
 
           integer:: ia,ja,ixyz,jxyz,ihl0,ihl1,ihl2,jj
           real(8):: xi(3),xj(3),xji(3),rij(3),rji(3),dji,sji(6),sii&
                ,hl2i,hl2j,tmp2i,tmp2j,hl1i,hl1j,tmp1i,tmp1j, fi(3)
 
           vatom(1:6,1:namax) = 0d0
           if( nl.eq.1 ) then
              do ia=1,natm
                 xi(1:3)= ra(1:3,ia)
                 do jj=1,lspr(0,ia)
                    ja= lspr(jj,ia)
                    xj(1:3)= ra(1:3,ja)
                    xji(1:3)= xj(1:3)-xi(1:3)
                    rji(1:3)= h(1:3,1)*xji(1) +h(1:3,2)*xji(2) +h(1:3,3)*xji(3)
                    rij(1:3)= -rji(1:3)
                    dji= sqrt(rji(1)**2 +rji(2)**2 +rji(3)**2)
                    if( dji.ge.rc ) cycle
                    do ihl1=1,nhl(1)
                       hl1i= hl1(ihl1,ia)
                       tmp1i= wgt12(ihl1)*hl1i*(1d0-hl1i)
                       do ihl0=1,nhl(0)
                          fi(1) = tmp1i*wgt11(ihl0,ihl1)*dgsf(1,ihl0,jj,ia)
                          fi(2) = tmp1i*wgt11(ihl0,ihl1)*dgsf(2,ihl0,jj,ia)
                          fi(3) = tmp1i*wgt11(ihl0,ihl1)*dgsf(3,ihl0,jj,ia)
                          sji(1) = -rji(1)* 0.5 * fi(1)
                          sji(2) = -rji(1)* 0.5 * fi(2)
                          sji(3) = -rji(1)* 0.5 * fi(3)
 
                          sji(4) = -0.25*(rji(1) * fi(2) + rji(2) * fi(1))
                          sji(5) = -0.25*(rji(1) * fi(3) + rji(3) * fi(1))
                          sji(6) = -0.25*(rji(2) * fi(3) + rji(3) * fi(2))
 
                          vatom(1,ia) = vatom(1,ia) + sji(1)
                          vatom(2,ia) = vatom(2,ia) + sji(2)
                          vatom(3,ia) = vatom(3,ia) + sji(3)
                          vatom(4,ia) = vatom(4,ia) + sji(4)
                          vatom(5,ia) = vatom(5,ia) + sji(5)
                          vatom(6,ia) = vatom(6,ia) + sji(6)
                          vatom(1,ja) = vatom(1,ja) + sji(1)
                          vatom(2,ja) = vatom(2,ja) + sji(2)
                          vatom(3,ja) = vatom(3,ja) + sji(3)
                          vatom(4,ja) = vatom(4,ja) + sji(4)
                          vatom(5,ja) = vatom(5,ja) + sji(5)
                          vatom(6,ja) = vatom(6,ja) + sji(6)
                       enddo
                    enddo
                 enddo
              enddo
           else if( nl.eq.2 ) then
              do ia=1,natm
                 xi(1:3)= ra(1:3,ia)
                 do jj=1,lspr(0,ia)
                    ja= lspr(jj,ia)
                    xj(1:3)= ra(1:3,ja)
                    xji(1:3)= xj(1:3)-xi(1:3)
                    rji(1:3)= h(1:3,1)*xji(1) +h(1:3,2)*xji(2) +h(1:3,3)*xji(3)
                    rij(1:3)= -rji(1:3)
                    dji= sqrt(rji(1)**2 +rji(2)**2 +rji(3)**2)
                    if( dji.ge.rc ) cycle
                    do ihl2=1,nhl(2)
                       hl2i= hl2(ihl2,ia)
                       tmp2i= wgt23(ihl2) *hl2i*(1d0-hl2i)
                       do ihl1=1,nhl(1)
                          hl1i= hl1(ihl1,ia)
                          tmp1i= wgt22(ihl1,ihl2) *hl1i*(1d0-hl1i)
                          do ihl0=1,nhl(0)
                             fi(1) = tmp1i*tmp2i*wgt21(ihl0,ihl1)*dgsf(1,ihl0,jj,ia)
                             fi(2) = tmp1i*tmp2i*wgt21(ihl0,ihl1)*dgsf(2,ihl0,jj,ia)
                             fi(3) = tmp1i*tmp2i*wgt21(ihl0,ihl1)*dgsf(3,ihl0,jj,ia)
 
                             sji(1) = -rji(1)* 0.5 * fi(1)
                             sji(2) = -rji(1)* 0.5 * fi(2)
                             sji(3) = -rji(1)* 0.5 * fi(3)
 
                             sji(4) = -0.25*(rji(1) * fi(2) + rji(2) * fi(1))
                             sji(5) = -0.25*(rji(1) * fi(3) + rji(3) * fi(1))
                             sji(6) = -0.25*(rji(2) * fi(3) + rji(3) * fi(2))
 
                             vatom(1,ia) = vatom(1,ia) + sji(1)
                             vatom(2,ia) = vatom(2,ia) + sji(2)
                             vatom(3,ia) = vatom(3,ia) + sji(3)
                             vatom(4,ia) = vatom(4,ia) + sji(4)
                             vatom(5,ia) = vatom(5,ia) + sji(5)
                             vatom(6,ia) = vatom(6,ia) + sji(6)
                             vatom(1,ja) = vatom(1,ja) + sji(1)
                             vatom(2,ja) = vatom(2,ja) + sji(2)
                             vatom(3,ja) = vatom(3,ja) + sji(3)
                             vatom(4,ja) = vatom(4,ja) + sji(4)
                             vatom(5,ja) = vatom(5,ja) + sji(5)
                             vatom(6,ja) = vatom(6,ja) + sji(6)
                          enddo
                       enddo
                    enddo
                 enddo
              enddo
           endif
 
      end subroutine compute_stress
 
 
 end module NN
 
 !-----------------------------------------------------------------------
 !     Local Variables:
 !     compile-command: "make pmd"
 !     End: