change for lipid last gly
[unres.git] / source / unres / src_MD-M / energy_p_new_barrier.F
index a20ec7b..443471d 100644 (file)
@@ -141,7 +141,7 @@ cmc Sep-06: egb takes care of dynamic ss bonds too
 cmc
 c      if (dyn_ss) call dyn_set_nss
 
-c      print *,"Processor",myrank," computed USCSC"
+C      print *,"Processor",myrank," computed USCSC"
 #ifdef TIMING
       time01=MPI_Wtime() 
 #endif
@@ -232,7 +232,7 @@ C#undef DEBUG
       enddo
 #endif
       endif
-c      print *,"Processor",myrank," left VEC_AND_DERIV"
+C      print *,"Processor",myrank," left VEC_AND_DERIV"
       if (ipot.lt.6) then
 #ifdef SPLITELE
          if (welec.gt.0d0.or.wvdwpp.gt.0d0.or.wel_loc.gt.0d0.or.
@@ -254,11 +254,11 @@ c      print *,"Processor",myrank," left VEC_AND_DERIV"
             eello_turn4=0.0d0
          endif
       else
-        write (iout,*) "Soft-spheer ELEC potential"
+C        write (iout,*) "Soft-spheer ELEC potential"
 c        call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3,
 c     &   eello_turn4)
       endif
-c      print *,"Processor",myrank," computed UELEC"
+C      print *,"Processor",myrank," computed UELEC"
 C
 C Calculate excluded-volume interaction energy between peptide groups
 C and side chains.
@@ -281,9 +281,9 @@ c
 C 
 C Calculate the disulfide-bridge and other energy and the contributions
 C from other distance constraints.
-cd    print *,'Calling EHPB'
+C      print *,'Calling EHPB'
       call edis(ehpb)
-cd    print *,'EHPB exitted succesfully.'
+C      print *,'EHPB exitted succesfully.'
 C
 C Calculate the virtual-bond-angle energy.
 C
@@ -300,13 +300,13 @@ C energy function
         ebe=0
         ethetacnstr=0
       endif
-c      print *,"Processor",myrank," computed UB"
+C       print *,"Processor",myrank," computed UB"
 C
 C Calculate the SC local energy.
 C
 C      print *,"TU DOCHODZE?"
       call esc(escloc)
-c      print *,"Processor",myrank," computed USC"
+C      print *,"Processor",myrank," computed USC"
 C
 C Calculate the virtual-bond torsional energy.
 C
@@ -325,7 +325,7 @@ C energy function
        etors=0
        edihcnstr=0
       endif
-c      print *,"Processor",myrank," computed Utor"
+C      print *,"Processor",myrank," computed Utor"
 C
 C 6/23/01 Calculate double-torsional energy
 C
@@ -334,7 +334,7 @@ C
       else
        etors_d=0
       endif
-c      print *,"Processor",myrank," computed Utord"
+C      print *,"Processor",myrank," computed Utord"
 C
 C 21/5/07 Calculate local sicdechain correlation energy
 C
@@ -344,7 +344,7 @@ C
         esccor=0.0d0
       endif
 C      print *,"PRZED MULIt"
-c      print *,"Processor",myrank," computed Usccorr"
+C      print *,"Processor",myrank," computed Usccorr"
 C 
 C 12/1/95 Multi-body terms
 C
@@ -3692,13 +3692,13 @@ c        if ((zmedi.gt.((0.5d0)*boxzsize)).or.
 c     &       (zmedi.lt.((-0.5d0)*boxzsize))) then
 c        go to 196
 c        endif
-          xmedi=mod(xmedi,boxxsize)
+          xmedi=dmod(xmedi,boxxsize)
           if (xmedi.lt.0) xmedi=xmedi+boxxsize
-          ymedi=mod(ymedi,boxysize)
+          ymedi=dmod(ymedi,boxysize)
           if (ymedi.lt.0) ymedi=ymedi+boxysize
-          zmedi=mod(zmedi,boxzsize)
+          zmedi=dmod(zmedi,boxzsize)
           if (zmedi.lt.0) zmedi=zmedi+boxzsize
-          zmedi2=mod(zmedi,boxzsize)
+          zmedi2=dmod(zmedi,boxzsize)
           if (zmedi2.lt.0) zmedi2=zmedi2+boxzsize
        if ((zmedi2.gt.bordlipbot)
      &.and.(zmedi2.lt.bordliptop)) then
@@ -3757,11 +3757,11 @@ c     &  .or. itype(i-1).eq.ntyp1
         xmedi=c(1,i)+0.5d0*dxi
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*dzi
-          xmedi=mod(xmedi,boxxsize)
+          xmedi=dmod(xmedi,boxxsize)
           if (xmedi.lt.0) xmedi=xmedi+boxxsize
-          ymedi=mod(ymedi,boxysize)
+          ymedi=dmod(ymedi,boxysize)
           if (ymedi.lt.0) ymedi=ymedi+boxysize
-          zmedi=mod(zmedi,boxzsize)
+          zmedi=dmod(zmedi,boxzsize)
           if (zmedi.lt.0) zmedi=zmedi+boxzsize
        if ((zmedi.gt.bordlipbot)
      &.and.(zmedi.lt.bordliptop)) then
@@ -3965,6 +3965,7 @@ C lipbufthick is thickenes of lipid buffore
        enddo
        enddo
        if (isubchap.eq.1) then
+C          print *,i,j
           xj=xj_temp-xmedi
           yj=yj_temp-ymedi
           zj=zj_temp-zmedi
@@ -4173,10 +4174,14 @@ C           print *,"bafter", gelc_long(1,i), gelc_long(1,j)
 C Lipidic part for lipscale
             gelc_long(3,j)=gelc_long(3,j)+
      &     ssgradlipj*eesij/2.0d0*lipscale**2
-
+C           if ((ssgradlipj*eesij/2.0d0*lipscale**2).ne.0.0 )
+C     &     write(iout,*) "WTF",j
             gelc_long(3,i)=gelc_long(3,i)+
      &     ssgradlipi*eesij/2.0d0*lipscale**2
 
+C            if ((ssgradlipi*eesij/2.0d0*lipscale**2).ne.0.0 )
+C     &     write(iout,*) "WTF",i
+
 *
 * Loop over residues i+1 thru j-1.
 *
@@ -4647,8 +4652,8 @@ c     &   a33*gmuji2(4)
 #endif
 cd          write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
 
-          if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
-     &            'eelloc',i,j,eel_loc_ij
+          if (energy_dec) write (iout,'(a6,2i5,0pf7.3,2f7.3)')
+     &            'eelloc',i,j,eel_loc_ij,a22*muij(1),a23*muij(2)
 c           if (eel_loc_ij.ne.0)
 c     &      write (iout,'(a4,2i4,8f9.5)')'chuj',
 c     &     i,j,a22,muij(1),a23,muij(2),a32,muij(3),a33,muij(4)
@@ -5120,6 +5125,7 @@ C Derivatives in gamma(i+1)
      &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
 
 C Cartesian derivatives
+!DIR$ UNROLL(0)
         do l=1,3
 c            ghalf1=0.5d0*agg(l,1)
 c            ghalf2=0.5d0*agg(l,2)
@@ -5809,6 +5815,7 @@ C
       include 'COMMON.CONTROL'
       include 'COMMON.SPLITELE'
       dimension ggg(3)
+      integer xshift,yshift,zshift
       evdw2=0.0D0
       evdw2_14=0.0d0
 c        print *,boxxsize,boxysize,boxzsize,'wymiary pudla'
@@ -6390,8 +6397,9 @@ C
 c      time11=dexp(-2*time)
 c      time12=1.0d0
       etheta=0.0D0
-c     write (*,'(a,i2)') 'EBEND ICG=',icg
+     write (*,'(a,i2)') 'EBEND ICG=',icg
       do i=ithet_start,ithet_end
+        if (i.le.2) cycle
         if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
      &  .or.itype(i).eq.ntyp1) cycle
 C Zero the energy function and its derivative at 0 or pi.
@@ -6410,7 +6418,8 @@ C Zero the energy function and its derivative at 0 or pi.
           ichir22=isign(1,itype(i))
          endif
 
-        if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
+        if (i.gt.3 ) then
+         if (itype(i-3).ne.ntyp1) then
 #ifdef OSF
          phii=phi(i)
           if (phii.ne.phii) phii=150.0
@@ -6423,6 +6432,11 @@ C Zero the energy function and its derivative at 0 or pi.
           y(1)=0.0D0
           y(2)=0.0D0
         endif
+        else 
+          y(1)=0.0D0
+          y(2)=0.0D0
+        endif
+
         if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
 #ifdef OSF
          phii1=phi(i+1)
@@ -6664,6 +6678,7 @@ C
       logical lprn /.false./, lprn1 /.false./
       etheta=0.0D0
       do i=ithet_start,ithet_end
+               if (i.le.2) cycle
 c        print *,i,itype(i-1),itype(i),itype(i-2)
 C        if (itype(i-1).eq.ntyp1) cycle
         if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
@@ -6681,7 +6696,8 @@ C        print *,i,theta(i)
           sinkt(k)=dsin(k*theti2)
         enddo
 C        print *,ethetai
-        if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
+        if (i.gt.3) then 
+         if (itype(i-3).ne.ntyp1) then
 #ifdef OSF
           phii=phi(i)
           if (phii.ne.phii) phii=150.0
@@ -6702,6 +6718,15 @@ C propagation of chirality for glycine type
             sinph1(k)=0.0d0
           enddo 
         endif
+        else
+          phii=0.0d0
+          do k=1,nsingle
+          ityp1=ithetyp((itype(i-2)))
+            cosph1(k)=0.0d0
+            sinph1(k)=0.0d0
+          enddo
+        endif
+
         if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
 #ifdef OSF
           phii1=phi(i+1)
@@ -8223,6 +8248,9 @@ c   3 = SC...Ca...Ca...SCi
           esccor=esccor+v1ij*cosphi+v2ij*sinphi
           gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
         enddo
+        if (energy_dec) write(iout,'(a9,2i4,f8.3,3i4)') "esccor",i,j,
+     & esccor,intertyp,
+     & isccori, isccori1
 c      write (iout,*) "EBACK_SC_COR",i,v1ij*cosphi+v2ij*sinphi,intertyp
         gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci
         if (lprn)
@@ -11569,7 +11597,8 @@ C--bordlipbot--buffore ends
       eliptran=0.0
       do i=ilip_start,ilip_end
 C       do i=1,1
-        if (itype(i).eq.ntyp1) cycle
+        if ((itype(i).eq.ntyp1).or.(itype(i+1).eq.ntyp1).or.(i.eq.nres))
+     &    cycle
 
         positi=(mod(((c(3,i)+c(3,i+1))/2.0d0),boxzsize))
         if (positi.le.0.0) positi=positi+boxzsize
@@ -12390,7 +12419,7 @@ C THIS FRAGMENT MAKES TUBE FINITE
 C       print *,mod(c(3,i+nres),boxzsize),bordlipbot,bordliptop
 c for each residue check if it is in lipid or lipid water border area
 C       respos=mod(c(3,i+nres),boxzsize)
-       print *,positi,bordtubebot,buftubebot,bordtubetop
+C       print *,positi,bordtubebot,buftubebot,bordtubetop
        if ((positi.gt.bordtubebot)
      & .and.(positi.lt.bordtubetop)) then
 C the energy transfer exist
@@ -12400,7 +12429,7 @@ C the energy transfer exist
 C lipbufthick is thickenes of lipid buffore
          sstube=sscalelip(fracinbuf)
          ssgradtube=-sscagradlip(fracinbuf)/tubebufthick
-         print *,ssgradtube, sstube,tubetranene(itype(i))
+C         print *,ssgradtube, sstube,tubetranene(itype(i))
          enetube(i)=enetube(i)+sstube*tubetranenepep
 C         gg_tube_SC(3,i)=gg_tube_SC(3,i)
 C     &+ssgradtube*tubetranene(itype(i))
@@ -12478,7 +12507,8 @@ C THIS FRAGMENT MAKES TUBE FINITE
 C       print *,mod(c(3,i+nres),boxzsize),bordlipbot,bordliptop
 c for each residue check if it is in lipid or lipid water border area
 C       respos=mod(c(3,i+nres),boxzsize)
-       print *,positi,bordtubebot,buftubebot,bordtubetop
+C       print *,positi,bordtubebot,buftubebot,bordtubetop
+
        if ((positi.gt.bordtubebot)
      & .and.(positi.lt.bordtubetop)) then
 C the energy transfer exist
@@ -12488,7 +12518,7 @@ C the energy transfer exist
 C lipbufthick is thickenes of lipid buffore
          sstube=sscalelip(fracinbuf)
          ssgradtube=-sscagradlip(fracinbuf)/tubebufthick
-         print *,ssgradtube, sstube,tubetranene(itype(i))
+C         print *,ssgradtube, sstube,tubetranene(itype(i))
          enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i))
 C         gg_tube_SC(3,i)=gg_tube_SC(3,i)
 C     &+ssgradtube*tubetranene(itype(i))
@@ -12608,17 +12638,17 @@ C now calculate distance from center of tube and direction vectors
       zmin=boxzsize
 
         do j=-1,1
-         vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize)
+         vectube(1)=dmod((c(1,i)+c(1,i+1))/2.0d0,boxxsize)
          vectube(1)=vectube(1)+boxxsize*j
-         vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize)
+         vectube(2)=dmod((c(2,i)+c(2,i+1))/2.0d0,boxysize)
          vectube(2)=vectube(2)+boxysize*j
-         vectube(3)=mod((c(3,i)+c(3,i+1))/2.0d0,boxzsize)
+         vectube(3)=dmod((c(3,i)+c(3,i+1))/2.0d0,boxzsize)
          vectube(3)=vectube(3)+boxzsize*j
 
 
-         xminact=abs(vectube(1)-tubecenter(1))
-         yminact=abs(vectube(2)-tubecenter(2))
-         zminact=abs(vectube(3)-tubecenter(3))
+         xminact=dabs(vectube(1)-tubecenter(1))
+         yminact=dabs(vectube(2)-tubecenter(2))
+         zminact=dabs(vectube(3)-tubecenter(3))
 
            if (xmin.gt.xminact) then
             xmin=xminact
@@ -12671,13 +12701,13 @@ C go to 667
          enecavtube(i)=0.0
          faccav=0.0
          else
-         denominator=(1.0+dcavtubpep*rdiff6*rdiff6)
+         denominator=(1.0d0+dcavtubpep*rdiff6*rdiff6)
          enecavtube(i)=
-     &   (bcavtubpep*rdiff+acavtubpep*sqrt(rdiff)+ccavtubpep)
+     &   (bcavtubpep*rdiff+acavtubpep*dsqrt(rdiff)+ccavtubpep)
      &   /denominator
          enecavtube(i)=0.0
-         faccav=((bcavtubpep*1.0d0+acavtubpep/2.0d0/sqrt(rdiff))
-     &   *denominator-(bcavtubpep*rdiff+acavtubpep*sqrt(rdiff)
+         faccav=((bcavtubpep*1.0d0+acavtubpep/2.0d0/dsqrt(rdiff))
+     &   *denominator-(bcavtubpep*rdiff+acavtubpep*dsqrt(rdiff)
      &   +ccavtubpep)*rdiff6**2.0d0/rdiff*dcavtubpep*12.0d0)
      &   /denominator**2.0d0
 C         faccav=0.0
@@ -12698,7 +12728,7 @@ C now direction of gg_tube vector
         enddo
 
        do i=itube_start,itube_end
-        enecavtube(i)=0.0 
+        enecavtube(i)=0.0d0
 C Lets not jump over memory as we use many times iti
          iti=itype(i)
 C lets ommit dummy atoms for now
@@ -12710,17 +12740,17 @@ C      .or.(iti.eq.10)
       ymin=boxysize
       zmin=boxzsize
         do j=-1,1
-         vectube(1)=mod((c(1,i+nres)),boxxsize)
+         vectube(1)=dmod((c(1,i+nres)),boxxsize)
          vectube(1)=vectube(1)+boxxsize*j
-         vectube(2)=mod((c(2,i+nres)),boxysize)
+         vectube(2)=dmod((c(2,i+nres)),boxysize)
          vectube(2)=vectube(2)+boxysize*j
-         vectube(3)=mod((c(3,i+nres)),boxzsize)
+         vectube(3)=dmod((c(3,i+nres)),boxzsize)
          vectube(3)=vectube(3)+boxzsize*j
 
 
-         xminact=abs(vectube(1)-tubecenter(1))
-         yminact=abs(vectube(2)-tubecenter(2))
-         zminact=abs(vectube(3)-tubecenter(3))
+         xminact=dabs(vectube(1)-tubecenter(1))
+         yminact=dabs(vectube(2)-tubecenter(2))
+         zminact=dabs(vectube(3)-tubecenter(3))
 
            if (xmin.gt.xminact) then
             xmin=xminact
@@ -12769,16 +12799,16 @@ C now direction of gg_tube vector
 C Now cavity term E=a(x+bsqrt(x)+c)/(1+dx^12)
          if (acavtub(iti).eq.0.0d0) then
 C go to 667
-         enecavtube(i+nres)=0.0
-         faccav=0.0
+         enecavtube(i+nres)=0.0d0
+         faccav=0.0d0
          else
-         denominator=(1.0+dcavtub(iti)*rdiff6*rdiff6)
+         denominator=(1.0d0+dcavtub(iti)*rdiff6*rdiff6)
          enecavtube(i+nres)=
-     &   (bcavtub(iti)*rdiff+acavtub(iti)*sqrt(rdiff)+ccavtub(iti))
+     &   (bcavtub(iti)*rdiff+acavtub(iti)*dsqrt(rdiff)+ccavtub(iti))
      &   /denominator
 C         enecavtube(i)=0.0
-         faccav=((bcavtub(iti)*1.0d0+acavtub(iti)/2.0d0/sqrt(rdiff))
-     &   *denominator-(bcavtub(iti)*rdiff+acavtub(iti)*sqrt(rdiff)
+         faccav=((bcavtub(iti)*1.0d0+acavtub(iti)/2.0d0/dsqrt(rdiff))
+     &   *denominator-(bcavtub(iti)*rdiff+acavtub(iti)*dsqrt(rdiff)
      &   +ccavtub(iti))*rdiff6**2.0d0/rdiff*dcavtub(iti)*12.0d0)
      &   /denominator**2.0d0
 C         faccav=0.0