projects / unres.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
Merge branch 'AFM' into multichain
[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 0044ec5..3192819 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -273,6 +273,11 @@ C      print *,"przed lipidami"
         call Eliptransfer(eliptran)
       endif
 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
 #endif
@@ -315,6 +320,7 @@ C
       energia(20)=Uconst+Uconst_back
       energia(21)=esccor
       energia(22)=eliptran
+      energia(23)=Eafmforce
 c    Here are the energies showed per procesor if the are more processors 
 c    per molecule then we sum it up in sum_energy subroutine 
 c      print *," Processor",myrank," calls SUM_ENERGY"
@@ -407,13 +413,14 @@ cMS$ATTRIBUTES C ::  proc_proc
       Uconst=energia(20)
       esccor=energia(21)
       eliptran=energia(22)
+      Eafmforce=energia(23)
 #ifdef SPLITELE
       etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1
      & +wang*ebe+wtor*etors+wscloc*escloc
      & +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5
      & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
      & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
-     & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran
+     & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+Eafmforce
 #else
       etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1)
      & +wang*ebe+wtor*etors+wscloc*escloc
@@ -421,6 +428,7 @@ cMS$ATTRIBUTES C ::  proc_proc
      & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
      & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
      & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran
+     & +Eafmforce
 #endif
       energia(0)=etot
 c detecting NaNQ
@@ -524,6 +532,7 @@ c      enddo
      &                wturn6*gcorr6_turn_long(j,i)+
      &                wstrain*ghpbc(j,i)
      &                +wliptran*gliptranc(j,i)
+     &                +gradafm(j,i)
 
         enddo
       enddo 
@@ -541,6 +550,8 @@ c      enddo
      &                wturn6*gcorr6_turn_long(j,i)+
      &                wstrain*ghpbc(j,i)
      &                +wliptran*gliptranc(j,i)
+     &                +gradafm(j,i)
+
         enddo
       enddo 
 #endif
@@ -677,6 +688,7 @@ c      enddo
      &                wsccor*gsccorc(j,i)
      &               +wscloc*gscloc(j,i)
      &               +wliptran*gliptranc(j,i)
+     &                +gradafm(j,i)
 #else
           gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+
      &                wel_loc*gel_loc(j,i)+
@@ -697,6 +709,8 @@ c      enddo
      &                wsccor*gsccorc(j,i)
      &               +wscloc*gscloc(j,i)
      &               +wliptran*gliptranc(j,i)
+     &                +gradafm(j,i)
+
 #endif
           gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+
      &                  wbond*gradbx(j,i)+
@@ -993,6 +1007,7 @@ C------------------------------------------------------------------------
       Uconst=energia(20)
       esccor=energia(21)
       eliptran=energia(22)
+      Eafmforce=energia(23) 
 #ifdef SPLITELE
       write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,
      &  estr,wbond,ebe,wang,
@@ -1001,7 +1016,7 @@ C------------------------------------------------------------------------
      &  ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
      &  eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,
      &  edihcnstr,ebr*nss,
-     &  Uconst,eliptran,wliptran,etot
+     &  Uconst,eliptran,wliptran,Eafmforce,etot
    10 format (/'Virtual-chain energies:'//
      & 'EVDW=  ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
      & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
@@ -1026,7 +1041,9 @@ C------------------------------------------------------------------------
      & 'ESS=   ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
      & 'UCONST= ',1pE16.6,' (Constraint energy)'/ 
      & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/
+     & 'EAFM=  ',1pE16.6,' (atomic-force microscopy)'/
      & 'ETOT=  ',1pE16.6,' (total)')
+
 #else
       write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,
      &  estr,wbond,ebe,wang,
@@ -1034,7 +1051,7 @@ C------------------------------------------------------------------------
      &  ecorr,wcorr,
      &  ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
      &  eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr,
-     &  ebr*nss,Uconst,eliptran,wliptran,etot
+     &  ebr*nss,Uconst,eliptran,wliptran,Eafmforc,etot
    10 format (/'Virtual-chain energies:'//
      & 'EVDW=  ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
      & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
@@ -1058,6 +1075,7 @@ C------------------------------------------------------------------------
      & 'ESS=   ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
      & 'UCONST=',1pE16.6,' (Constraint energy)'/ 
      & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/
+     & 'EAFM=  ',1pE16.6,' (atomic-force microscopy)'/
      & 'ETOT=  ',1pE16.6,' (total)')
 #endif
       return
@@ -3267,11 +3285,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)
@@ -3295,6 +3320,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
@@ -3359,6 +3388,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
@@ -3413,6 +3446,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
@@ -10224,3 +10261,80 @@ C       eliptran=elpitran+0.0 ! I am in water
        enddo
        return
        end
+C---------------------------------------------------------
+C AFM soubroutine for constant force
+       subroutine AFMforce(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)
+      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=-forceAFMconst*(dist-distafminit)
+      do i=1,3
+      gradafm(i,afmend-1)=-forceAFMconst*diffafm(i)/dist
+      gradafm(i,afmbeg-1)=forceAFMconst*diffafm(i)/dist
+      enddo
+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
+
UNRESPACK 3.x