anchor.diff:      !$OMP PARALLEL DO PRIVATE(ni)
       do i = 1, nm2
           ni = n2(i)
  -        do m_1 = 1, ni
  +        do m_1 = 2, ni
               do i_1 = 1, ni
                   if (x2(i,i_1,1,m_1) < cut_distance) then
  -                    ksi2(i, i_1, m_1) = 1.0d0 - sin(x2(i,i_1,1,m_1) * inv_cut)
  +                    ! ksi2(i, i_1, m_1) = 1.0d0 - sin(x2(i,i_1,1,m_1) * inv_cut)
  +                    ksi2(i, i_1, m_1) = 1.0d0 / x2(i,i_1,1,m_1)**6 &
  +                    & * x2(i,i_1,2,1) * x2(i,i_1,2,m_1) 
                   endif
               enddo
           enddo
       enddo
       !$OMP END PARALLEL DO
  +
  +    allocate(ksi31(nm1, maxval(n1), maxval(n1), maxval(n1)))
  +    allocate(ksi32(nm2, maxval(n2), maxval(n2), maxval(n2)))
  +
  +    ksi31 = 0.0d0
  +    ksi32 = 0.0d0
  +
  +    
  +    do i = 1, nm1
  +        ni = n1(i)
  +        do m_2 = 2, ni
  +            do m_1 = 2, ni
  +                do i_1 = 1, ni
  +
  +                    if (m_2.eq.m_1) cycle 
  +
  +                    di = x1(i, i_1, 1, m_1)
  +                    dj = x1(i, i_1, 1, m_2)
  +                    angleij = x1(i, i_1, 3 + m_1, m_2) 
  +
  +                    dij = sqrt(di**2 + dj**2 - 2.0d0*di*dj*cos(angleij))
  +
  +                    !cosi = sqrt(1.0d0 - (dj * sin(angleij)/dij)**2) 
  +                    !cosj = sqrt(1.0d0 - (di * sin(angleij)/dij)**2) 
  +
  +                    cosi = (dj**2 + dij**2 - di**2)/(2.0d0 * dj * dij)
  +                    cosj = (di**2 + dij**2 - dj**2)/(2.0d0 * di * dij)
  +
  +                    ! ksi31(i,i_1,m_1,m_2) = (1.0d0 + cos(angleij) * cosi * cosj) / (di * dj * dij)**3
  +
  +                    ksi31(i,i_1,m_1,m_2) = x1(i,i_1,2,1) * x1(i,i_1,2,m_1) * x1(i,i_1,2,m_2) * &
  +                    & (1.0d0 + cos(angleij) * cosi * cosj) / (di * dj * dij)**3
  +
  +                    ! write(*,*) "KSI31", i, i_1, m_1, m_2, ksi31(i, i_1, m_1, m_2)
  +                enddo
  +            enddo
  +        enddo
  +    enddo
  +
  +
  +    do i = 1, nm2
  +        ni = n2(i)
  +        do m_2 = 2, ni
  +            do m_1 = 2, ni
  +                do i_1 = 1, ni
  +
  +                    if (m_2.eq.m_1) cycle 
  +
  +                    di = x2(i, i_1, 1, m_1)
  +                    dj = x2(i, i_1, 1, m_2)
  +                    angleij = x2(i, i_1, 3 + m_1, m_2) 
  +
  +                    dij = sqrt(di**2 + dj**2 - 2.0d0*di*dj*cos(angleij))
  +
  +                    ! cosi = sqrt(1.0d0 - (dj * sin(angleij)/dij)**2) 
  +                    ! cosj = sqrt(1.0d0 - (di * sin(angleij)/dij)**2) 
  +
  +                    cosi = (dj**2 + dij**2 - di**2)/(2.0d0 * dj * dij)
  +                    cosj = (di**2 + dij**2 - dj**2)/(2.0d0 * di * dij)
  +
  +                    ! ksi32(i,i_1,m_1,m_2) = (1.0d0 + cos(angleij) * cosi * cosj) / (di * dj * dij)**3
  +
  +                    ksi32(i,i_1,m_1,m_2) = x2(i,i_1,2,1) * x2(i,i_1,2,m_1) * x2(i,i_1,2,m_2) * &
  +                    & (1.0d0 + cos(angleij) * cosi * cosj) / (di * dj * dij)**3
  +
  +                enddo
  +            enddo
  +        enddo
  +    enddo
  +
   
       allocate(cosp1(nm1, maxval(n1), order, maxval(n1), pmax1))
       allocate(sinp1(nm1, maxval(n1), order, maxval(n1), pmax1))
  
positive.code:     !$OMP PARALLEL DO PRIVATE(ni)
      do i = 1, nm1
          ni = n1(i)
          do m_1 = 1, ni
              do i_1 = 1, ni
                  if (x1(i,i_1,1,m_1) < cut_distance) then
                      ksi1(i, i_1, m_1) = 1.0d0 - sin(x1(i,i_1,1,m_1) * inv_cut)
                  endif
              enddo
          enddo
  
negative.code: 
                  enddo
  
                  aadist = aadist + d * (ksi1(m_1) * ksi2(m_2) + angular &
                  & * angular_scale)
                  !aadist = aadist + d * (ksi1(m_1) * ksi2(m_2))
  
              end if
          end do
  