corrections to src_MD-M for the same energy as src_MD
[unres.git] / source / unres / src_MD-M / energy_p_new_barrier.F
index b9d17fe..9ea982a 100644 (file)
@@ -85,7 +85,7 @@ C FG slaves receive the WEIGHTS array
         time_Bcastw=time_Bcastw+MPI_Wtime()-time00
 c        call chainbuild_cart
       endif
-c      print *,'Processor',myrank,' calling etotal ipot=',ipot
+C      print *,'Processor',myrank,' calling etotal ipot=',ipot
 c      print *,'Processor',myrank,' nnt=',nnt,' nct=',nct
 #else
 c      if (modecalc.eq.12.or.modecalc.eq.14) then
@@ -98,6 +98,7 @@ c      endif
 C 
 C Compute the side-chain and electrostatic interaction energy
 C
+C      write(iout,*) "zaczynam liczyc energie"
       goto (101,102,103,104,105,106) ipot
 C Lennard-Jones potential.
   101 call elj(evdw)
@@ -117,10 +118,13 @@ C Gay-Berne-Vorobjev potential (shifted LJ, angular dependence).
       goto 107
 C Soft-sphere potential
   106 call e_softsphere(evdw)
+C      write(iout,*) "skonczylem ipoty"
+
 C
 C Calculate electrostatic (H-bonding) energy of the main chain.
 C
   107 continue
+C           write(iout,*) "skonczylem ipoty"
 cmc
 cmc Sep-06: egb takes care of dynamic ss bonds too
 cmc
@@ -442,7 +446,6 @@ cMS$ATTRIBUTES C ::  proc_proc
 #endif
       double precision gradbufc(3,maxres),gradbufx(3,maxres),
      &  glocbuf(4*maxres),gradbufc_sum(3,maxres),gloc_scbuf(3,maxres)
-#endif
       include 'COMMON.SETUP'
       include 'COMMON.IOUNITS'
       include 'COMMON.FFIELD'
@@ -715,7 +718,7 @@ c      enddo
         do i=1,4*nres
           glocbuf(i)=gloc(i,icg)
         enddo
-#define DEBUG
+#undef DEBUG
 #ifdef DEBUG
       write (iout,*) "gloc_sc before reduce"
       do i=1,nres
@@ -744,7 +747,7 @@ c      enddo
         call MPI_Reduce(gloc_scbuf(1,1),gloc_sc(1,1,icg),3*nres,
      &    MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
         time_reduce=time_reduce+MPI_Wtime()-time00
-#define DEBUG
+#undef DEBUG
 #ifdef DEBUG
       write (iout,*) "gloc_sc after reduce"
       do i=1,nres
@@ -2338,6 +2341,7 @@ C
         endif
 c        if (i.gt. iatel_s+2 .and. i.lt.iatel_e+5) then
         if (i.gt. nnt+2 .and. i.lt.nct+2) then
+c        write(iout,*) (itype(i-2))
           iti = itortyp(itype(i-2))
         else
           iti=ntortyp+1
@@ -2718,7 +2722,7 @@ C
       dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
      &          erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
       double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
-     &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
+     &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),eel_loc_ij
       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
      &    dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
      &    num_conti,j1,j2
@@ -2761,6 +2765,7 @@ c        call vec_and_deriv
         time01=MPI_Wtime()
 #endif
         call set_matrices
+c        write (iout,*) "after set matrices"
 #ifdef TIMING
         time_mat=time_mat+MPI_Wtime()-time01
 #endif
@@ -2797,6 +2802,7 @@ c 9/27/08 AL Split the interaction loop to ensure load balancing of turn terms
 C
 C Loop over i,i+2 and i,i+3 pairs of the peptide groups
 C
+c      write(iout,*) "przed turnem3 loop"
       do i=iturn3_start,iturn3_end
         if (itype(i).eq.21 .or. itype(i+1).eq.21 
      &  .or. itype(i+2).eq.21 .or. itype(i+3).eq.21) cycle
@@ -2889,7 +2895,7 @@ C-------------------------------------------------------------------------------
       dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
      &          erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
       double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
-     &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
+     &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),a22,a23,a32,a33
       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
      &    dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
      &    num_conti,j1,j2
@@ -3659,7 +3665,7 @@ c        write (iout,*) "eturn4 i",i," j",j," j1",j1," j2",j2
         iti1=itortyp(itype(i+1))
         iti2=itortyp(itype(i+2))
         iti3=itortyp(itype(i+3))
-c        write(iout,*) "iti1",iti1," iti2",iti2," iti3",iti3
+C        write(iout,*) i,"iti1",iti1," iti2",iti2," iti3",iti3,itype(i+3)
         call transpose2(EUg(1,1,i+1),e1t(1,1))
         call transpose2(Eug(1,1,i+2),e2t(1,1))
         call transpose2(Eug(1,1,i+3),e3t(1,1))
@@ -4540,7 +4546,8 @@ C
       logical lprn /.false./, lprn1 /.false./
       etheta=0.0D0
       do i=ithet_start,ithet_end
-        if (itype(i-1).eq.21) cycle
+        if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or.
+     &(itype(i).eq.ntyp1)) cycle
         dethetai=0.0d0
         dephii=0.0d0
         dephii1=0.0d0
@@ -4550,7 +4557,8 @@ C
           coskt(k)=dcos(k*theti2)
           sinkt(k)=dsin(k*theti2)
         enddo
-        if (i.gt.3 .and. itype(i-2).ne.21) then
+C        if (i.gt.3) then
+        if (i.gt.3 .and. itype(imax0(i-3,1)).ne.ntyp1) then
 #ifdef OSF
           phii=phi(i)
           if (phii.ne.phii) phii=150.0
           enddo
         else
           phii=0.0d0
-          ityp1=nthetyp+1
+          ityp1=ithetyp(itype(i-2))
           do k=1,nsingle
             cosph1(k)=0.0d0
             sinph1(k)=0.0d0
           enddo 
         endif
-        if (i.lt.nres .and. itype(i).ne.21) then
+        if ((i.lt.nres).and. itype(i+1).ne.ntyp1) then
 #ifdef OSF
           phii1=phi(i+1)
           if (phii1.ne.phii1) phii1=150.0
@@ -4585,7 +4593,7 @@ C
           enddo
         else
           phii1=0.0d0
-          ityp3=nthetyp+1
+          ityp3=ithetyp(itype(i))
           do k=1,nsingle
             cosph2(k)=0.0d0
             sinph2(k)=0.0d0
@@ -4695,6 +4703,8 @@ C
      &   i,theta(i)*rad2deg,phii*rad2deg,
      &   phii1*rad2deg,ethetai
         etheta=etheta+ethetai
+        if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
+     &      'ebend',i,ethetai
         if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii
         if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1
         gloc(nphi+i-2,icg)=wang*dethetai
@@ -5553,7 +5563,8 @@ c     lprn=.true.
       etors=0.0D0
       do i=iphi_start,iphi_end
         if (itype(i-2).eq.21 .or. itype(i-1).eq.21 
-     &       .or. itype(i).eq.21) cycle
+     &       .or. itype(i).eq.21
+     &       .or. itype(i-3).eq.ntyp1) cycle
       etors_ii=0.0D0
         itori=itortyp(itype(i-2))
         itori1=itortyp(itype(i-1))
@@ -5641,15 +5652,18 @@ C 6/23/01 Compute double torsional energy
       include 'COMMON.IOUNITS'
       include 'COMMON.FFIELD'
       include 'COMMON.TORCNSTR'
+      include 'COMMON.CONTROL'
       logical lprn
 C Set lprn=.true. for debugging
       lprn=.false.
 c     lprn=.true.
       etors_d=0.0D0
-c      write(iout,*) "a tu??"
+C      write(iout,*) "a tu??"
       do i=iphid_start,iphid_end
         if (itype(i-2).eq.21 .or. itype(i-1).eq.21
-     &      .or. itype(i).eq.21 .or. itype(i+1).eq.21) cycle
+     &      .or. itype(i).eq.21 .or. itype(i+1).eq.21
+     &       .or. itype(i-3).eq.ntyp1) cycle
+        etors_d_ii=0.0D0
         itori=itortyp(itype(i-2))
         itori1=itortyp(itype(i-1))
         itori2=itortyp(itype(i))
@@ -5669,6 +5683,8 @@ C Regular cosine and sine terms
           sinphi2=dsin(j*phii1)
           etors_d=etors_d+v1cij*cosphi1+v1sij*sinphi1+
      &     v2cij*cosphi2+v2sij*sinphi2
+          if (energy_dec) etors_d_ii=etors_d_ii+
+     &     v1cij*cosphi1+v1sij*sinphi1+v2cij*cosphi2+v2sij*sinphi2
           gloci1=gloci1+j*(v1sij*cosphi1-v1cij*sinphi1)
           gloci2=gloci2+j*(v2sij*cosphi2-v2cij*sinphi2)
         enddo
@@ -5684,12 +5700,17 @@ C Regular cosine and sine terms
             sinphi1m2=dsin(l*phii-(k-l)*phii1)
             etors_d=etors_d+v1cdij*cosphi1p2+v2cdij*cosphi1m2+
      &        v1sdij*sinphi1p2+v2sdij*sinphi1m2
+            if (energy_dec) etors_d_ii=etors_d_ii+
+     &        v1cdij*cosphi1p2+v2cdij*cosphi1m2+
+     &        v1sdij*sinphi1p2+v2sdij*sinphi1m2
             gloci1=gloci1+l*(v1sdij*cosphi1p2+v2sdij*cosphi1m2
      &        -v1cdij*sinphi1p2-v2cdij*sinphi1m2)
             gloci2=gloci2+(k-l)*(v1sdij*cosphi1p2-v2sdij*cosphi1m2
      &        -v1cdij*sinphi1p2+v2cdij*sinphi1m2) 
           enddo
         enddo
+          if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
+     &         'etor_d',i,etors_d_ii
         gloc(i-3,icg)=gloc(i-3,icg)+wtor_d*gloci1
         gloc(i-2,icg)=gloc(i-2,icg)+wtor_d*gloci2
       enddo