emomo correction
[unres4.git] / source / unres / energy.F90
index 4ae8d1b..fce83dd 100644 (file)
        integer ishield_listbuf(-1:nres), &
        shield_listbuf(maxcontsshi,-1:nres),k,j,i,iii,impishi,mojint,jjj
 !       print *,"I START ENERGY"
-       imatupdate=1
+       imatupdate=100
 !       if (mod(itime_mat,imatupdate).eq.0) call make_SCSC_inter_list
 !      real(kind=8),  dimension(:),allocatable::  fac_shieldbuf 
 !      real(kind=8), dimension(:,:,:),allocatable:: &
         time_Bcastw=time_Bcastw+MPI_Wtime()-time00
 !        call chainbuild_cart
       endif
-!       if (mod(itime_mat,imatupdate).eq.0) call make_SCSC_inter_list
+!       print *,"itime_mat",itime_mat,imatupdate
+        if (nfgtasks.gt.1) then 
+        call MPI_Bcast(itime_mat,1,MPI_INT,king,FG_COMM,IERROR)
+        endif
+       if (mod(itime_mat,imatupdate).eq.0) call make_SCp_inter_list
+       if (mod(itime_mat,imatupdate).eq.0) call make_SCSC_inter_list
+       if (mod(itime_mat,imatupdate).eq.0) call make_pp_inter_list
 
 !      print *,'Processor',myrank,' calling etotal ipot=',ipot
 !      print *,'Processor',myrank,' nnt=',nnt,' nct=',nct
        etors_nucl=0.0d0
        estr_nucl=0.0d0
        ecorr3_nucl=0.0d0
+       ecorr_nucl=0.0d0
        ebe_nucl=0.0d0
        evdwsb=0.0d0
        eelsb=0.0d0
        eelpsb=0.0d0
        evdwpp=0.0d0
        eespp=0.0d0
+       etors_d_nucl=0.0d0
       endif
 !      write(iout,*) ecorr_nucl,"ecorr_nucl",nres_molec(2)
-      print *,"before ecatcat",wcatcat
+!      print *,"before ecatcat",wcatcat
+      if (nres_molec(5).gt.0) then
       if (nfgtasks.gt.1) then
       if (fg_rank.eq.0) then
       call ecatcat(ecationcation)
       else
       call ecats_prot_amber(ecation_prot)
       endif
-      if (nres_molec(2).gt.0) then
+      else
+      ecationcation=0.0d0
+      ecation_prot=0.0d0
+      endif
+      if ((nres_molec(2).gt.0).and.(nres_molec(1).gt.0)) then
       call eprot_sc_base(escbase)
       call epep_sc_base(epepbase)
       call eprot_sc_phosphate(escpho)
 !      include 'COMMON.SBRIDGE'
       logical :: lprn
 !el local variables
-      integer :: iint,itypi,itypi1,itypj,subchap
+      integer :: iint,itypi,itypi1,itypj,subchap,icont
       real(kind=8) :: rrij,xi,yi,zi,sig,rij_shift,fac,e1,e2,sigm,epsi
       real(kind=8) :: evdw,sig0ij
       real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
       dPOLdOM1=0.0d0
 
 
-      do i=iatsc_s,iatsc_e
+      do icont=g_listscsc_start,g_listscsc_end
+      i=newcontlisti(icont)
+      j=newcontlistj(icont)
+
+!      do i=iatsc_s,iatsc_e
 !C        print *,"I am in EVDW",i
         itypi=iabs(itype(i,1))
 !        if (i.ne.47) cycle
 !
 ! Calculate SC interaction energy.
 !
-        do iint=1,nint_gr(i)
-          do j=istart(i,iint),iend(i,iint)
+!        do iint=1,nint_gr(i)
+!          do j=istart(i,iint),iend(i,iint)
             IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
               call dyn_ssbond_ene(i,j,evdwij)
               evdw=evdw+evdwij
 ! Calculate angular part of the gradient.
             call sc_grad
             ENDIF    ! dyn_ss            
-          enddo      ! j
-        enddo        ! iint
+!          enddo      ! j
+!        enddo        ! iint
       enddo          ! i
 !       print *,"ZALAMKA", evdw
 !      write (iout,*) "Number of loop steps in EGB:",ind
 #endif
 #else
         if (i.gt. nnt+2 .and. i.lt.nct+2) then
-!         write(iout,*) "i,",molnum(i)
+!         write(iout,*) "i,",molnum(i),nloctyp
 !         print *, "i,",molnum(i),i,itype(i-2,1)
         if (molnum(i).eq.1) then
+          if (itype(i-2,1).eq.ntyp1) then
+           iti=nloctyp
+          else
           iti = itype2loc(itype(i-2,1))
+          endif
         else
           iti=nloctyp
         endif
                                              0.0d0,1.0d0,0.0d0,&
                                              0.0d0,0.0d0,1.0d0/),shape(unmat)) 
 !el local variables
-      integer :: i,k,j
+      integer :: i,k,j,icont
       real(kind=8) :: ees,evdw1,eel_loc,eello_turn3,eello_turn4
       real(kind=8) :: fac,t_eelecij,fracinbuf
     
 ! Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
 !
 !      print *,"iatel_s,iatel_e,",iatel_s,iatel_e
-      do i=iatel_s,iatel_e
+!      do i=iatel_s,iatel_e
+! JPRDLC
+       do icont=g_listpp_start,g_listpp_end
+        i=newcontlistppi(icont)
+        j=newcontlistppj(icont)
         if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
         dxi=dc(1,i)
         dyi=dc(2,i)
 
 !        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
         num_conti=num_cont_hb(i)
-        do j=ielstart(i),ielend(i)
+!        do j=ielstart(i),ielend(i)
 !          write (iout,*) i,j,itype(i,1),itype(j,1)
           if (itype(j,1).eq.ntyp1.or. itype(j+1,1).eq.ntyp1) cycle
           call eelecij(i,j,ees,evdw1,eel_loc)
-        enddo ! j
+!        enddo ! j
         num_cont_hb(i)=num_conti
       enddo   ! i
 !      write (iout,*) "Number of loop steps in EELEC:",ind
 !             sss_ele_grad=0.0d0
 !            print *,sss_ele_cut,sss_ele_grad,&
 !            (rij),r_cut_ele,rlamb_ele
-!            if (sss_ele_cut.le.0.0) go to 128
+            if (sss_ele_cut.le.0.0) go to 128
 
           rmij=1.0D0/rij
           r3ij=rrmij*rmij
 !grad            enddo
 !grad          enddo
 ! 9/28/08 AL Gradient compotents will be summed only at the end
-          ggg(1)=facvdw*xj &
+          ggg(1)=facvdw*xj+sss_ele_grad*rmij*evdwij*xj &
            *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
-          ggg(2)=facvdw*yj &
+          ggg(2)=facvdw*yj+sss_ele_grad*rmij*evdwij*yj &
            *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
-          ggg(3)=facvdw*zj &
+          ggg(3)=facvdw*zj+sss_ele_grad*rmij*evdwij*zj &
            *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
 
           do k=1,3
           +a32*gmuij1(3)&
           +a33*gmuij1(4))&
          *fac_shield(i)*fac_shield(j)&
-                    *sss_ele_cut
+                    *sss_ele_cut     &
+         *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
 
 !c         write(iout,*) "derivative over thatai"
 !c         write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
          gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+&
            geel_loc_ij*wel_loc&
          *fac_shield(i)*fac_shield(j)&
-                    *sss_ele_cut
+                    *sss_ele_cut &
+         *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
 
 
 !c  Derivative over j residue
         gloc(nphi+j,icg)=gloc(nphi+j,icg)+&
            geel_loc_ji*wel_loc&
          *fac_shield(i)*fac_shield(j)&
-                    *sss_ele_cut
+                    *sss_ele_cut &
+         *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
 
 
          geel_loc_ji=&
          gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+&
            geel_loc_ji*wel_loc&
          *fac_shield(i)*fac_shield(j)&
-                    *sss_ele_cut
+                    *sss_ele_cut &
+         *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
 #endif
 
 !          write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
                 ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij) &
                      *sss_ele_cut &
                      *fac_shield(i)*fac_shield(j)
+!                     *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
 
                 ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij) &
                      *sss_ele_cut &
                      *fac_shield(i)*fac_shield(j)
+!                     *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
 
 ! Diagnostics. Comment out or remove after debugging!
 !               ees0p(num_conti,i)=0.5D0*fac3*ees0pij
                   gacontp_hb1(k,num_conti,i)= & !ghalfp+
                     (ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i)) &
                    + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) &
-                     *sss_ele_cut*fac_shield(i)*fac_shield(j)
+                     *sss_ele_cut*fac_shield(i)*fac_shield(j) ! &
+!                     *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
 
                   gacontp_hb2(k,num_conti,i)= & !ghalfp+
                     (ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j)) &
                    + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)&
-                     *sss_ele_cut*fac_shield(i)*fac_shield(j)
+                     *sss_ele_cut*fac_shield(i)*fac_shield(j)!   &
+!                     *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
 
                   gacontp_hb3(k,num_conti,i)=gggp(k) &
                      *sss_ele_cut*fac_shield(i)*fac_shield(j)
+!                     *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
 
                   gacontm_hb1(k,num_conti,i)= & !ghalfm+
                     (ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i)) &
                    + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) &
                      *sss_ele_cut*fac_shield(i)*fac_shield(j)
+!                     *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
 
                   gacontm_hb2(k,num_conti,i)= & !ghalfm+
                     (ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j)) &
                    + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) &
                      *sss_ele_cut*fac_shield(i)*fac_shield(j)
+!                     *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
 
                   gacontm_hb3(k,num_conti,i)=gggm(k) &
                      *sss_ele_cut*fac_shield(i)*fac_shield(j)
+!                     *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
 
                 enddo
 ! Diagnostics. Comment out or remove after debugging!
 !C Derivatives in theta
         gloc(nphi+i,icg)=gloc(nphi+i,icg) &
        +0.5d0*(gpizda1(1,1)+gpizda1(2,2))*wturn3&
-        *fac_shield(i)*fac_shield(j)
+        *fac_shield(i)*fac_shield(j) &
+        *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
         gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)&
        +0.5d0*(gpizda2(1,1)+gpizda2(2,2))*wturn3&
-        *fac_shield(i)*fac_shield(j)
+        *fac_shield(i)*fac_shield(j) &
+        *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
+
 !C#endif
 
 
 !      include 'COMMON.CONTROL'
       real(kind=8),dimension(3) :: ggg
 !el local variables
-      integer :: i,iint,j,k,iteli,itypj,subchap
+      integer :: i,iint,j,k,iteli,itypj,subchap,icont
       real(kind=8) :: evdw2,evdw2_14,xi,yi,zi,xj,yj,zj,rrij,fac,&
                    e1,e2,evdwij,rij
       real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
       evdw2_14=0.0d0
 !d    print '(a)','Enter ESCP'
 !d    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
-      do i=iatscp_s,iatscp_e
+!      do i=iatscp_s,iatscp_e
+       do icont=g_listscp_start,g_listscp_end
+        i=newcontlistscpi(icont)
+        j=newcontlistscpj(icont)
         if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
         iteli=itel(i)
         xi=0.5D0*(c(1,i)+c(1,i+1))
           zi=mod(zi,boxzsize)
           if (zi.lt.0) zi=zi+boxzsize
 
-        do iint=1,nscp_gr(i)
+!        do iint=1,nscp_gr(i)
 
-        do j=iscpstart(i,iint),iscpend(i,iint)
+!        do j=iscpstart(i,iint),iscpend(i,iint)
           itypj=iabs(itype(j,1))
           if (itypj.eq.ntyp1) cycle
 ! Uncomment following three lines for SC-p interactions
             gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
             gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
           enddo
-        enddo
+!        enddo
 
-        enddo ! iint
+!        enddo ! iint
       enddo ! i
       do i=1,nct
         do j=1,3
 !     &   dscp1,dscp2,sumene
 !        sumene = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
         escloc = escloc + sumene
+       if (energy_dec) write (2,*) "i",i," itype",itype(i,1)," it",it, &
+        " escloc",sumene,escloc,it,itype(i,1)
 !        write (2,*) "i",i," escloc",sumene,escloc,it,itype(i,1)
 !     & ,zz,xx,yy
 !#define DEBUG
       allocate(uygrad(3,3,2,nres))
       allocate(uzgrad(3,3,2,nres))
 !(3,3,2,maxres)
+! allocateion of lists JPRDLA
+      allocate(newcontlistppi(200*nres))
+      allocate(newcontlistscpi(200*nres))
+      allocate(newcontlisti(200*nres))
+      allocate(newcontlistppj(200*nres))
+      allocate(newcontlistscpj(200*nres))
+      allocate(newcontlistj(200*nres))
 
       return
       end subroutine alloc_ener_arrays
 !        r06 = rcat0**6
 !        r012 = r06**2
 !        k0 = 332*(2*2)/80
-        Evan1cat=epscalc*(r012/rcal**6)
-        Evan2cat=epscalc*2*(r06/rcal**3)
+        Evan1cat=epscalc*(r012/(rcal**6))
+        Evan2cat=epscalc*2*(r06/(rcal**3))
         Eeleccat=k0/ract
         r7 = rcal**7
         r4 = rcal**4
           gradcatcat(k,i)=gradcatcat(k,i)-gg(k)
           gradcatcat(k,j)=gradcatcat(k,j)+gg(k)
         enddo
-
+        if (energy_dec) write (iout,*) i,j,Evan1cat,Evan2cat,Eeleccat,&
+         r012,rcal**6,ichargecat(itypi)*ichargecat(itypj)
 !        write(iout,*) "ecatcat",i,j, ecationcation,xj,yj,zj
         ecationcation=ecationcation+Evan1cat+Evan2cat+Eeleccat
        enddo
 !c! No charges - do nothing
            eheadtail = 0.0d0
 
-          ELSE IF (isel.eq.1 .and. iabs(Qj).eq.1) THEN
+          ELSE IF (isel.eq.1) THEN
 !c! Nonpolar-charge interactions
           if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
             Qi=Qi*2
            eheadtail = epol
 !           eheadtail = 0.0d0
 
-          ELSE IF (isel.eq.3 .and. icharge(itypj).eq.2) THEN
+          ELSE IF (isel.eq.3) THEN
 !c! Dipole-charge interactions
           if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
             Qi=Qi*2
             Qj=Qj*2
             Qij=Qij*2
            endif
+           write(iout,*) "KURWA0",d1
+
            CALL edq_cat(ecl, elj, epol)
           eheadtail = ECL + elj + epol
 !           eheadtail = 0.0d0
 
-          ELSE IF ((isel.eq.2.and.   &
-               iabs(Qi).eq.1).and.  &
-               nstatecat(itypi,itypj).eq.1) THEN
+          ELSE IF ((isel.eq.2)) THEN
 
 !c! Same charge-charge interaction ( +/+ or -/- )
           if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
        real(kind=8) ::  facd4, adler, Fgb, facd3
        integer troll,jj,istate
        real (kind=8) :: dcosom1(3),dcosom2(3)
+       evdw=0.0d0
        eps_out=80.0d0
        sss_ele_cut=1.0d0
 !       print *,"EVDW KURW",evdw,nres
 
 !c!-------------------------------------------------------------------
 !c! ecl
+       write(iout,*) "KURWA2",Rhead
        sparrow  = w1 * Qj * om1
        hawk     = w2 * Qj * Qj * (1.0d0 - sqom2)
        ECL = sparrow / Rhead**2.0d0 &
       include 'mpif.h'
       real*8 :: xi,yi,zi,xj,yj,zj,xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp
       real*8 :: dist_init, dist_temp,r_buff_list
-      integer:: contlisti(50*nres),contlistj(50*nres)
-      integer :: newcontlisti(50*nres),newcontlistj(50*nres) 
+      integer:: contlisti(200*nres),contlistj(200*nres)
+!      integer :: newcontlisti(200*nres),newcontlistj(200*nres) 
       integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint,ilist_sc,g_ilist_sc
       integer displ(0:nprocs),i_ilist_sc(0:nprocs),ierr
 !            print *,"START make_SC"
+          r_buff_list=5.0
             ilist_sc=0
             do i=iatsc_s,iatsc_e
              itypi=iabs(itype(i,1))
       write (iout,*) i,newcontlisti(i),newcontlistj(i)
       enddo
 #endif
-
+        call int_bounds(g_ilist_sc,g_listscsc_start,g_listscsc_end)
       return
       end subroutine make_SCSC_inter_list
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
       subroutine make_SCp_inter_list
+      use MD_data,  only: itime_mat
+
       include 'mpif.h'
       real*8 :: xi,yi,zi,xj,yj,zj,xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp
       real*8 :: dist_init, dist_temp,r_buff_list
-      integer:: contlistscpi(50*nres),contlistscpj(50*nres)
-      integer :: newcontlistscpi(50*nres),newcontlistscpj(50*nres)
+      integer:: contlistscpi(200*nres),contlistscpj(200*nres)
+!      integer :: newcontlistscpi(200*nres),newcontlistscpj(200*nres)
       integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint,ilist_scp,g_ilist_scp
       integer displ(0:nprocs),i_ilist_scp(0:nprocs),ierr
-            print *,"START make_SC"
+!            print *,"START make_SC"
+      r_buff_list=5.0
             ilist_scp=0
       do i=iatscp_s,iatscp_e
         if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
           yj=yj_safe-yi
           zj=zj_safe-zi
        endif
+#ifdef DEBUG
+                ! r_buff_list is a read value for a buffer 
+               if ((sqrt(dist_init).le.(r_cut_ele)).and.(ifirstrun.eq.0)) then
+! Here the list is created
+                 ilist_scp_first=ilist_scp_first+1
+! this can be substituted by cantor and anti-cantor
+                 contlistscpi_f(ilist_scp_first)=i
+                 contlistscpj_f(ilist_scp_first)=j
+              endif
+#endif
 ! r_buff_list is a read value for a buffer 
                if (sqrt(dist_init).le.(r_cut_ele+r_buff_list)) then
 ! Here the list is created
       do i=1,g_ilist_scp
       write (iout,*) i,newcontlistscpi(i),newcontlistscpj(i)
       enddo
+
+!      if (ifirstrun.eq.0) ifirstrun=1
+!      do i=1,ilist_scp_first
+!       do j=1,g_ilist_scp
+!        if ((newcontlistscpi(j).eq.contlistscpi_f(i)).and.&
+!         (newcontlistscpj(j).eq.contlistscpj_f(i))) go to 126
+!        enddo
+!       print *,itime_mat,"ERROR matrix needs updating"
+!       print *,contlistscpi_f(i),contlistscpj_f(i)
+!  126  continue
+!      enddo
 #endif
+        call int_bounds(g_ilist_scp,g_listscp_start,g_listscp_end)
+
       return
       end subroutine make_SCp_inter_list
 
 !-----------------------------------------------------------------------------
 
 
+      subroutine make_pp_inter_list
+      include 'mpif.h'
+      real*8 :: xi,yi,zi,xj,yj,zj,xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp
+      real*8 :: xmedj,ymedj,zmedj
+      real*8 :: dist_init, dist_temp,r_buff_list,dxi,dyi,dzi,xmedi,ymedi,zmedi
+      real*8 :: dx_normi,dy_normi,dz_normi,dxj,dyj,dzj,dx_normj,dy_normj,dz_normj
+      integer:: contlistppi(200*nres),contlistppj(200*nres)
+!      integer :: newcontlistppi(200*nres),newcontlistppj(200*nres)
+      integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint,ilist_pp,g_ilist_pp
+      integer displ(0:nprocs),i_ilist_pp(0:nprocs),ierr
+!            print *,"START make_SC"
+            ilist_pp=0
+      r_buff_list=5.0
+      do i=iatel_s,iatel_e
+        if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
+        dxi=dc(1,i)
+        dyi=dc(2,i)
+        dzi=dc(3,i)
+        dx_normi=dc_norm(1,i)
+        dy_normi=dc_norm(2,i)
+        dz_normi=dc_norm(3,i)
+        xmedi=c(1,i)+0.5d0*dxi
+        ymedi=c(2,i)+0.5d0*dyi
+        zmedi=c(3,i)+0.5d0*dzi
+          xmedi=dmod(xmedi,boxxsize)
+          if (xmedi.lt.0) xmedi=xmedi+boxxsize
+          ymedi=dmod(ymedi,boxysize)
+          if (ymedi.lt.0) ymedi=ymedi+boxysize
+          zmedi=dmod(zmedi,boxzsize)
+          if (zmedi.lt.0) zmedi=zmedi+boxzsize
+             do j=ielstart(i),ielend(i)
+!          write (iout,*) i,j,itype(i,1),itype(j,1)
+          if (itype(j,1).eq.ntyp1.or. itype(j+1,1).eq.ntyp1) cycle
+! 1,j)
+          dxj=dc(1,j)
+          dyj=dc(2,j)
+          dzj=dc(3,j)
+          dx_normj=dc_norm(1,j)
+          dy_normj=dc_norm(2,j)
+          dz_normj=dc_norm(3,j)
+!          xj=c(1,j)+0.5D0*dxj-xmedi
+!          yj=c(2,j)+0.5D0*dyj-ymedi
+!          zj=c(3,j)+0.5D0*dzj-zmedi
+          xj=c(1,j)+0.5D0*dxj
+          yj=c(2,j)+0.5D0*dyj
+          zj=c(3,j)+0.5D0*dzj
+          xj=mod(xj,boxxsize)
+          if (xj.lt.0) xj=xj+boxxsize
+          yj=mod(yj,boxysize)
+          if (yj.lt.0) yj=yj+boxysize
+          zj=mod(zj,boxzsize)
+          if (zj.lt.0) zj=zj+boxzsize
+
+      dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
+      xj_safe=xj
+      yj_safe=yj
+      zj_safe=zj
+      do xshift=-1,1
+      do yshift=-1,1
+      do zshift=-1,1
+          xj=xj_safe+xshift*boxxsize
+          yj=yj_safe+yshift*boxysize
+          zj=zj_safe+zshift*boxzsize
+          dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
+          if(dist_temp.lt.dist_init) then
+            dist_init=dist_temp
+            xj_temp=xj
+            yj_temp=yj
+            zj_temp=zj
+          endif
+       enddo
+       enddo
+       enddo
+
+      if (sqrt(dist_init).le.(r_cut_ele+r_buff_list)) then
+! Here the list is created
+                 ilist_pp=ilist_pp+1
+! this can be substituted by cantor and anti-cantor
+                 contlistppi(ilist_pp)=i
+                 contlistppj(ilist_pp)=j
+              endif
+             enddo
+             enddo
+!             enddo
+#ifdef DEBUG
+      write (iout,*) "before MPIREDUCE",ilist_pp
+      do i=1,ilist_pp
+      write (iout,*) i,contlistppi(i),contlistppj(i)
+      enddo
+#endif
+      if (nfgtasks.gt.1)then
+
+        call MPI_Reduce(ilist_pp,g_ilist_pp,1,&
+          MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
+!        write(iout,*) "before bcast",g_ilist_sc
+        call MPI_Gather(ilist_pp,1,MPI_INTEGER,&
+                        i_ilist_pp,1,MPI_INTEGER,king,FG_COMM,IERR)
+        displ(0)=0
+        do i=1,nfgtasks-1,1
+          displ(i)=i_ilist_pp(i-1)+displ(i-1)
+        enddo
+!        write(iout,*) "before gather",displ(0),displ(1)
+        call MPI_Gatherv(contlistppi,ilist_pp,MPI_INTEGER,&
+                         newcontlistppi,i_ilist_pp,displ,MPI_INTEGER,&
+                         king,FG_COMM,IERR)
+        call MPI_Gatherv(contlistppj,ilist_pp,MPI_INTEGER,&
+                         newcontlistppj,i_ilist_pp,displ,MPI_INTEGER,&
+                         king,FG_COMM,IERR)
+        call MPI_Bcast(g_ilist_pp,1,MPI_INT,king,FG_COMM,IERR)
+!        write(iout,*) "before bcast",g_ilist_sc
+!        call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+        call MPI_Bcast(newcontlistppi,g_ilist_pp,MPI_INT,king,FG_COMM,IERR)
+        call MPI_Bcast(newcontlistppj,g_ilist_pp,MPI_INT,king,FG_COMM,IERR)
+
+!        call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+
+        else
+        g_ilist_pp=ilist_pp
+
+        do i=1,ilist_pp
+        newcontlistppi(i)=contlistppi(i)
+        newcontlistppj(i)=contlistppj(i)
+        enddo
+        endif
+        call int_bounds(g_ilist_pp,g_listpp_start,g_listpp_end)
+#ifdef DEBUG
+      write (iout,*) "after MPIREDUCE",g_ilist_pp
+      do i=1,g_ilist_pp
+      write (iout,*) i,newcontlistppi(i),newcontlistppj(i)
+      enddo
+#endif
+      return
+      end subroutine make_pp_inter_list
+
+!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
+
+
+
       end module energy