first changes for shielding

[unres.git] / source / unres / src_MD-M / energy_p_new_barrier.F

diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F
index 2bd6f22..9d6982d 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -137,6 +137,13 @@ c      print *,"Processor",myrank," computed USCSC"
 #ifdef TIMING
       time_vec=time_vec+MPI_Wtime()-time01
 #endif
+C Introduction of shielding effect first for each peptide group
+C the shielding factor is set this factor is describing how each
+C peptide group is shielded by side-chains
+C the matrix - shield_fac(i) the i index describe the ith between i and i+1
+      if (shield_mode.gt.0) then
+       call set_shield_fac
+      endif
 c      print *,"Processor",myrank," left VEC_AND_DERIV"
       if (ipot.lt.6) then
 #ifdef SPLITELE
@@ -275,6 +282,8 @@ C      print *,"przed lipidami"
 C      print *,"za lipidami"
       if (AFMlog.gt.0) then
         call AFMforce(Eafmforce)
+      else if (selfguide.gt.0) then
+        call AFMvel(Eafmforce)
       endif
 #ifdef TIMING
       time_enecalc=time_enecalc+MPI_Wtime()-time00
@@ -3283,11 +3292,18 @@ C 14/01/2014 TURN3,TUNR4 does no go under periodic boundry condition
         if (i.le.1) cycle
 C        write(iout,*) "tu jest i",i
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
+C changes suggested by Ana to avoid out of bounds
+     & .or.((i+4).gt.nres)
+     & .or.((i-1).le.0)
+C end of changes by Ana
      &  .or. itype(i+2).eq.ntyp1
-     &  .or. itype(i+3).eq.ntyp1
-     &  .or. itype(i-1).eq.ntyp1
-     &  .or. itype(i+4).eq.ntyp1
-     &  ) cycle
+     &  .or. itype(i+3).eq.ntyp1) cycle
+        if(i.gt.1)then
+          if(itype(i-1).eq.ntyp1)cycle
+        end if
+        if(i.LT.nres-3)then
+          if (itype(i+4).eq.ntyp1) cycle
+        end if
         dxi=dc(1,i)
         dyi=dc(2,i)
         dzi=dc(3,i)
@@ -3311,6 +3327,10 @@ C        write(iout,*) "tu jest i",i
       do i=iturn4_start,iturn4_end
         if (i.le.1) cycle
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
+C changes suggested by Ana to avoid out of bounds
+     & .or.((i+5).gt.nres)
+     & .or.((i-1).le.0)
+C end of changes suggested by Ana
      &    .or. itype(i+3).eq.ntyp1
      &    .or. itype(i+4).eq.ntyp1
      &    .or. itype(i+5).eq.ntyp1
@@ -3375,6 +3395,10 @@ c
       do i=iatel_s,iatel_e
         if (i.le.1) cycle
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
+C changes suggested by Ana to avoid out of bounds
+     & .or.((i+2).gt.nres)
+     & .or.((i-1).le.0)
+C end of changes by Ana
      &  .or. itype(i+2).eq.ntyp1
      &  .or. itype(i-1).eq.ntyp1
      &                ) cycle
@@ -3429,6 +3453,10 @@ c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
 C          write (iout,*) i,j
          if (j.le.1) cycle
           if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1
+C changes suggested by Ana to avoid out of bounds
+     & .or.((j+2).gt.nres)
+     & .or.((j-1).le.0)
+C end of changes by Ana
      & .or.itype(j+2).eq.ntyp1
      & .or.itype(j-1).eq.ntyp1
      &) cycle
@@ -10273,3 +10301,47 @@ C AFM soubroutine for constant force
 C      print *,'AFM',Eafmforce
       return
       end
+C---------------------------------------------------------
+C AFM subroutine with pseudoconstant velocity
+       subroutine AFMvel(Eafmforce)
+       implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.GEO'
+      include 'COMMON.VAR'
+      include 'COMMON.LOCAL'
+      include 'COMMON.CHAIN'
+      include 'COMMON.DERIV'
+      include 'COMMON.NAMES'
+      include 'COMMON.INTERACT'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CALC'
+      include 'COMMON.CONTROL'
+      include 'COMMON.SPLITELE'
+      include 'COMMON.SBRIDGE'
+      real*8 diffafm(3)
+C Only for check grad COMMENT if not used for checkgrad
+C      totT=3.0d0
+C--------------------------------------------------------
+C      print *,"wchodze"
+      dist=0.0d0
+      Eafmforce=0.0d0
+      do i=1,3
+      diffafm(i)=c(i,afmend)-c(i,afmbeg)
+      dist=dist+diffafm(i)**2
+      enddo
+      dist=dsqrt(dist)
+      Eafmforce=0.5d0*forceAFMconst
+     & *(distafminit+totTafm*velAFMconst-dist)**2
+C      Eafmforce=-forceAFMconst*(dist-distafminit)
+      do i=1,3
+      gradafm(i,afmend-1)=-forceAFMconst*
+     &(distafminit+totTafm*velAFMconst-dist)
+     &*diffafm(i)/dist
+      gradafm(i,afmbeg-1)=forceAFMconst*
+     &(distafminit+totTafm*velAFMconst-dist)
+     &*diffafm(i)/dist
+      enddo
+C      print *,'AFM',Eafmforce,totTafm*velAFMconst,dist
+      return
+      end
+