HOMOL klapaucjusz correction
[unres.git] / source / unres / src_MD / energy_p_new_barrier.F
index 115cb57..cbf930c 100644 (file)
@@ -3895,6 +3895,7 @@ C Derivatives in gamma(i+1)
         gel_loc_turn3(i+1)=gel_loc_turn3(i+1)
      &    +0.5d0*(pizda(1,1)+pizda(2,2))
 C Cartesian derivatives
+!DIR$ UNROLL(0)
         do l=1,3
 c            ghalf1=0.5d0*agg(l,1)
 c            ghalf2=0.5d0*agg(l,2)
@@ -4361,6 +4362,7 @@ C
       include 'COMMON.VAR'
       include 'COMMON.INTERACT'
       include 'COMMON.IOUNITS'
+      include 'COMMON.CONTROL'
       dimension ggg(3)
       ehpb=0.0D0
 cd      write(iout,*)'edis: nhpb=',nhpb,' fbr=',fbr
@@ -4396,6 +4398,16 @@ cd          write (iout,*) "eij",eij
         else if (ii.gt.nres .and. jj.gt.nres) then
 c Restraints from contact prediction
           dd=dist(ii,jj)
+          if (constr_dist.eq.11) then
+            ehpb=ehpb+fordepth(i)**4.0d0
+     &          *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i))
+            fac=fordepth(i)**4.0d0
+     &          *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd
+          if (energy_dec) write (iout,'(a6,2i5,f15.6,2f8.3)') 
+     &     "edisl",ii,jj,
+     &     fordepth(i)**4.0d0*rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)),
+     &     fordepth(i),dd
+           else
           if (dhpb1(i).gt.0.0d0) then
             ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
             fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd
@@ -4413,7 +4425,8 @@ C
 C Evaluate gradient.
 C
             fac=waga*rdis/dd
-          endif  
+          endif
+          endif
           do j=1,3
             ggg(j)=fac*(c(j,jj)-c(j,ii))
           enddo
@@ -4429,6 +4442,18 @@ C
 C Calculate the distance between the two points and its difference from the
 C target distance.
           dd=dist(ii,jj)
+         if (constr_dist.eq.11) then
+            ehpb=ehpb+fordepth(i)**4.0d0
+     &           *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i))
+            fac=fordepth(i)**4.0d0
+     &           *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd
+          if (energy_dec) write (iout,'(a6,2i5,f15.6,2f8.3)') 
+     7     "edisl",ii,jj,
+     &     fordepth(i)**4.0d0*rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)),
+     &     fordepth(i),dd
+c          if (energy_dec)
+c     &      write (iout,*) fac
+         else   
           if (dhpb1(i).gt.0.0d0) then
             ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
             fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd
@@ -4446,6 +4471,7 @@ C Evaluate gradient.
 C
             fac=waga*rdis/dd
           endif
+         endif
 cd      print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd,
 cd   &   ' waga=',waga,' fac=',fac
             do j=1,3
@@ -4471,7 +4497,10 @@ cgrad        enddo
           enddo
         endif
       enddo
-      ehpb=0.5D0*ehpb
+      if (constr_dist.ne.11) ehpb=0.5D0*ehpb
+c      do i=1,nres
+c        write (iout,*) "ghpbc",i,(ghpbc(j,i),j=1,3)
+c      enddo
       return
       end
 C--------------------------------------------------------------------------
@@ -4903,7 +4932,7 @@ c     &     " ithet_end",ithet_end
           sinkt(k)=dsin(k*theti2)
         enddo
 C        if (i.gt.3) then
-        if (i.gt.3 .and. itype(imax0(i-3,1)).ne.ntyp1) then
+        if (i.gt.3 .and. itype(max0(i-3,1)).ne.ntyp1) then
 #ifdef OSF
           phii=phi(i)
           if (phii.ne.phii) phii=150.0
@@ -6048,9 +6077,13 @@ c     write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
          j = jres_homo(ii)
          dij=dist(i,j)
 c        write (iout,*) "dij(",i,j,") =",dij
+         nexl=0
          do k=1,constr_homology
 c           write(iout,*) ii,k,i,j,l_homo(k,ii),dij,odl(k,ii)
-           if(.not.l_homo(k,ii)) cycle
+           if(.not.l_homo(k,ii)) then
+             nexl=nexl+1
+             cycle
+           endif
            distance(k)=odl(k,ii)-dij
 c          write (iout,*) "distance(",k,") =",distance(k)
 c
@@ -6069,6 +6102,7 @@ c
      &                     (distance(k)**2+sigma_odlir(k,ii)**2))
            endif
          enddo
+c         write (iout,*) "distance: ii",ii," nexl",nexl
          
 
 c         min_odl=minval(distancek)
@@ -6089,7 +6123,15 @@ c        write (iout,* )"min_odl",min_odl
          write (iout,*) "distancek",(distancek(k),k=1,constr_homology)
          write (iout,* )"min_odl",min_odl
 #endif
+#ifdef OLDRESTR
          odleg2=0.0d0
+#else
+         if (waga_dist.ge.0.0d0) then
+           odleg2=nexl
+         else 
+           odleg2=0.0d0
+         endif 
+#endif
          do k=1,constr_homology
 c Nie wiem po co to liczycie jeszcze raz!
 c            odleg3=-waga_dist(iset)*((distance(i,j,k)**2)/ 
@@ -6232,10 +6274,10 @@ c      write (iout,*) idihconstr_start_homo,idihconstr_end_homo
       enddo
 #endif
       do i=idihconstr_start_homo,idihconstr_end_homo
-        kat2=0.0d0
 c        betai=beta(i,i+1,i+2,i+3)
         betai = phi(i)
 c       write (iout,*) "betai =",betai
+        kat2=0.0d0
         do k=1,constr_homology
           dih_diff(k)=pinorm(dih(k,i)-betai)
 c         write (iout,*) "dih_diff(",k,") =",dih_diff(k)
@@ -6243,13 +6285,16 @@ c          if (dih_diff(i,k).gt.3.14159) dih_diff(i,k)=
 c     &                                   -(6.28318-dih_diff(i,k))
 c          if (dih_diff(i,k).lt.-3.14159) dih_diff(i,k)=
 c     &                                   6.28318+dih_diff(i,k)
-
+#ifdef OLD_DIHED
           kat3=-0.5d0*dih_diff(k)**2*sigma_dih(k,i) ! waga_angle rmvd from Gaussian argument
+#else
+          kat3=(dcos(dih_diff(k))-1)*sigma_dih(k,i) ! waga_angle rmvd from Gaussian argument
+#endif
 c         kat3=-0.5d0*waga_angle*dih_diff(k)**2*sigma_dih(k,i)
           gdih(k)=dexp(kat3)
           kat2=kat2+gdih(k)
-c          write(iout,*) "kat2=", kat2, "exp(kat3)=", exp(kat3)
-c          write(*,*)""
+c          write(iout,*) "i",i," k",k," sigma",sigma_dih(k,i),
+c     &     " kat2=", kat2, "gdih=",gdih(k)
         enddo
 c       write (iout,*) "gdih",(gdih(k),k=1,constr_homology) ! exps
 c       write (iout,*) "i",i," betai",betai," kat2",kat2 ! sum of exps
@@ -6271,7 +6316,11 @@ c ----------------------------------------------------------------------
         sum_gdih=kat2
         sum_sgdih=0.0d0
         do k=1,constr_homology
+#ifdef OLD_DIHED
           sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i)  ! waga_angle rmvd
+#else
+          sgdih=-gdih(k)*dsin(dih_diff(k))*sigma_dih(k,i)  ! waga_angle rmvd
+#endif
 c         sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i)*waga_angle
           sum_sgdih=sum_sgdih+sgdih
         enddo
@@ -6281,7 +6330,7 @@ c       grad_dih3=sum_sgdih/sum_gdih
 c      write(iout,*)i,k,gdih,sgdih,beta(i+1,i+2,i+3,i+4),grad_dih3
 ccc      write(iout,747) "GRAD_KAT_1", i, nphi, icg, grad_dih3,
 ccc     & gloc(nphi+i-3,icg)
-        gloc(i,icg)=gloc(i,icg)+grad_dih3
+        gloc(i-3,icg)=gloc(i-3,icg)+grad_dih3
 c        if (i.eq.25) then
 c        write(iout,*) "i",i,"icg",icg,"gloc(",i,icg,")",gloc(i,icg)
 c        endif
@@ -6350,7 +6399,9 @@ c
           utheta_i=-0.5d0*theta_diff(k)**2*sigma_theta(k,i) ! waga_theta rmvd from Gaussian argument
 c         utheta_i=-0.5d0*waga_theta*theta_diff(k)**2*sigma_theta(k,i) ! waga_theta?
           gtheta(k)=dexp(utheta_i) ! + min_utheta_i?
-          gutheta_i=gutheta_i+dexp(utheta_i)   ! Sum of Gaussians (pk)
+          gutheta_i=gutheta_i+gtheta(k)  ! Sum of Gaussians (pk)
+c          write (iout,*) "i",i," k",k," sigma_theta",sigma_theta(k,i),
+c     &     " gtheta",gtheta(k)
 c         Gradient for single Gaussian restraint in subr Econstr_back
 c         dutheta(j-2)=dutheta(j-2)+wfrag_back(1,i,iset)*dtheta_i/(ii-1)
 c
@@ -6422,7 +6473,9 @@ c
 c         usc_diff(i)=-0.5d0*waga_d*(dxx**2+dyy**2+dzz**2)*sigma_d(k,i) ! waga_d?
 c         uscdiffk(k)=usc_diff(i)
           guscdiff2(k)=dexp(usc_diff_i) ! without min_scdiff
-          guscdiff(i)=guscdiff(i)+dexp(usc_diff_i)   !Sum of Gaussians (pk)
+c          write(iout,*) "i",i," k",k," sigma_d",sigma_d(k,i),
+c     &       " guscdiff2",guscdiff2(k)
+          guscdiff(i)=guscdiff(i)+guscdiff2(k)  !Sum of Gaussians (pk)
 c          write (iout,'(i5,6f10.5)') j,xxtab(j),yytab(j),zztab(j),
 c     &      xxref(j),yyref(j),zzref(j)
         enddo