Merge branch 'devel' into AFM
[unres.git] / source / unres / src_Eshel / SRC-SURPLUS / energy_split-sep.F
diff --git a/source/unres/src_Eshel/SRC-SURPLUS/energy_split-sep.F b/source/unres/src_Eshel/SRC-SURPLUS/energy_split-sep.F
new file mode 100644 (file)
index 0000000..81e4d81
--- /dev/null
@@ -0,0 +1,476 @@
+      subroutine etotal_long(energia)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+c
+c Compute the long-range slow-varying contributions to the energy
+c
+#ifndef ISNAN
+      external proc_proc
+#ifdef WINPGI
+cMS$ATTRIBUTES C ::  proc_proc
+#endif
+#endif
+#ifdef MPI
+      include "mpif.h"
+      double precision weights_(n_ene)
+#endif
+      include 'COMMON.SETUP'
+      include 'COMMON.IOUNITS'
+      double precision energia(0:n_ene)
+      include 'COMMON.FFIELD'
+      include 'COMMON.DERIV'
+      include 'COMMON.INTERACT'
+      include 'COMMON.SBRIDGE'
+      include 'COMMON.CHAIN'
+      include 'COMMON.VAR'
+      include 'COMMON.LOCAL'
+      include 'COMMON.MD'
+c      write(iout,'(a,i2)')'Calling etotal_long ipot=',ipot
+      if (modecalc.eq.12.or.modecalc.eq.14) then
+#ifdef MPI
+c        if (fg_rank.eq.0) call int_from_cart1(.false.)
+#else
+        call int_from_cart1(.false.)
+#endif
+      endif
+#ifdef MPI      
+c      write(iout,*) "ETOTAL_LONG Processor",fg_rank,
+c     & " absolute rank",myrank," nfgtasks",nfgtasks
+      call flush(iout)
+      if (nfgtasks.gt.1) then
+        time00=MPI_Wtime()
+C FG slaves call the following matching MPI_Bcast in ERGASTULUM
+        if (fg_rank.eq.0) then
+          call MPI_Bcast(3,1,MPI_INTEGER,king,FG_COMM,IERROR)
+c          write (iout,*) "Processor",myrank," BROADCAST iorder"
+c          call flush(iout)
+C FG master sets up the WEIGHTS_ array which will be broadcast to the 
+C FG slaves as WEIGHTS array.
+          weights_(1)=wsc
+          weights_(2)=wscp
+          weights_(3)=welec
+          weights_(4)=wcorr
+          weights_(5)=wcorr5
+          weights_(6)=wcorr6
+          weights_(7)=wel_loc
+          weights_(8)=wturn3
+          weights_(9)=wturn4
+          weights_(10)=wturn6
+          weights_(11)=wang
+          weights_(12)=wscloc
+          weights_(13)=wtor
+          weights_(14)=wtor_d
+          weights_(15)=wstrain
+          weights_(16)=wvdwpp
+          weights_(17)=wbond
+          weights_(18)=scal14
+          weights_(21)=wsccor
+C FG Master broadcasts the WEIGHTS_ array
+          call MPI_Bcast(weights_(1),n_ene,
+     &        MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
+        else
+C FG slaves receive the WEIGHTS array
+          call MPI_Bcast(weights(1),n_ene,
+     &        MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
+          wsc=weights(1)
+          wscp=weights(2)
+          welec=weights(3)
+          wcorr=weights(4)
+          wcorr5=weights(5)
+          wcorr6=weights(6)
+          wel_loc=weights(7)
+          wturn3=weights(8)
+          wturn4=weights(9)
+          wturn6=weights(10)
+          wang=weights(11)
+          wscloc=weights(12)
+          wtor=weights(13)
+          wtor_d=weights(14)
+          wstrain=weights(15)
+          wvdwpp=weights(16)
+          wbond=weights(17)
+          scal14=weights(18)
+          wsccor=weights(21)
+        endif
+        call MPI_Bcast(dc(1,1),6*nres,MPI_DOUBLE_PRECISION,
+     &    king,FG_COMM,IERR)
+         time_Bcast=time_Bcast+MPI_Wtime()-time00
+         time_Bcastw=time_Bcastw+MPI_Wtime()-time00
+c        call chainbuild_cart
+c        call int_from_cart1(.false.)
+      endif
+c      write (iout,*) 'Processor',myrank,
+c     &  ' calling etotal_short ipot=',ipot
+c      call flush(iout)
+c      print *,'Processor',myrank,' nnt=',nnt,' nct=',nct
+#endif     
+cd    print *,'nnt=',nnt,' nct=',nct
+C
+C Compute the side-chain and electrostatic interaction energy
+C
+      goto (101,102,103,104,105,106) ipot
+C Lennard-Jones potential.
+  101 call elj_long(evdw)
+cd    print '(a)','Exit ELJ'
+      goto 107
+C Lennard-Jones-Kihara potential (shifted).
+  102 call eljk_long(evdw)
+      goto 107
+C Berne-Pechukas potential (dilated LJ, angular dependence).
+  103 call ebp_long(evdw)
+      goto 107
+C Gay-Berne potential (shifted LJ, angular dependence).
+  104 call egb_long(evdw,evdw_p,evdw_m)
+      goto 107
+C Gay-Berne-Vorobjev potential (shifted LJ, angular dependence).
+  105 call egbv_long(evdw)
+      goto 107
+C Soft-sphere potential
+  106 call e_softsphere(evdw)
+C
+C Calculate electrostatic (H-bonding) energy of the main chain.
+C
+  107 continue
+      call vec_and_deriv
+      if (ipot.lt.6) then
+#ifdef SPLITELE
+         if (welec.gt.0d0.or.wvdwpp.gt.0d0.or.wel_loc.gt.0d0.or.
+     &       wturn3.gt.0d0.or.wturn4.gt.0d0 .or. wcorr.gt.0.0d0
+     &       .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.d0
+     &       .or. wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0 ) then
+#else
+         if (welec.gt.0d0.or.wel_loc.gt.0d0.or.
+     &       wturn3.gt.0d0.or.wturn4.gt.0d0 .or. wcorr.gt.0.0d0
+     &       .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.d0
+     &       .or. wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0 ) then
+#endif
+           call eelec_scale(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
+         else
+            ees=0
+            evdw1=0
+            eel_loc=0
+            eello_turn3=0
+            eello_turn4=0
+         endif
+      else
+c        write (iout,*) "Soft-spheer ELEC potential"
+        call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3,
+     &   eello_turn4)
+      endif
+C
+C Calculate excluded-volume interaction energy between peptide groups
+C and side chains.
+C
+      if (ipot.lt.6) then
+       if(wscp.gt.0d0) then
+        call escp_long(evdw2,evdw2_14)
+       else
+        evdw2=0
+        evdw2_14=0
+       endif
+      else
+        call escp_soft_sphere(evdw2,evdw2_14)
+      endif
+C 
+C 12/1/95 Multi-body terms
+C
+      n_corr=0
+      n_corr1=0
+      if ((wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0 
+     &    .or. wturn6.gt.0.0d0) .and. ipot.lt.6) then
+         call multibody_eello(ecorr,ecorr5,ecorr6,eturn6,n_corr,n_corr1)
+c         write (2,*) 'n_corr=',n_corr,' n_corr1=',n_corr1,
+c     &" ecorr",ecorr," ecorr5",ecorr5," ecorr6",ecorr6," eturn6",eturn6
+      else
+         ecorr=0.0d0
+         ecorr5=0.0d0
+         ecorr6=0.0d0
+         eturn6=0.0d0
+      endif
+      if ((wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) .and. ipot.lt.6) then
+         call multibody_hb(ecorr,ecorr5,ecorr6,n_corr,n_corr1)
+      endif
+C 
+C If performing constraint dynamics, call the constraint energy
+C  after the equilibration time
+      if(usampl.and.totT.gt.eq_time) then
+         call EconstrQ   
+         call Econstr_back
+      else
+         Uconst=0.0d0
+         Uconst_back=0.0d0
+      endif
+C 
+C Sum the energies
+C
+      do i=1,n_ene
+        energia(i)=0.0d0
+      enddo
+      energia(1)=evdw
+#ifdef SCP14
+      energia(2)=evdw2-evdw2_14
+      energia(18)=evdw2_14
+#else
+      energia(2)=evdw2
+      energia(18)=0.0d0
+#endif
+#ifdef SPLITELE
+      energia(3)=ees
+      energia(16)=evdw1
+#else
+      energia(3)=ees+evdw1
+      energia(16)=0.0d0
+#endif
+      energia(4)=ecorr
+      energia(5)=ecorr5
+      energia(6)=ecorr6
+      energia(7)=eel_loc
+      energia(8)=eello_turn3
+      energia(9)=eello_turn4
+      energia(10)=eturn6
+      energia(20)=Uconst+Uconst_back
+      energia(22)=evdw_p
+      energia(23)=evdw_m
+      call sum_energy(energia,.true.)
+c      write (iout,*) "Exit ETOTAL_LONG"
+      call flush(iout)
+      return
+      end
+c------------------------------------------------------------------------------
+      subroutine etotal_short(energia)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+c
+c Compute the short-range fast-varying contributions to the energy
+c
+#ifndef ISNAN
+      external proc_proc
+#ifdef WINPGI
+cMS$ATTRIBUTES C ::  proc_proc
+#endif
+#endif
+#ifdef MPI
+      include "mpif.h"
+      double precision weights_(n_ene)
+#endif
+      include 'COMMON.SETUP'
+      include 'COMMON.IOUNITS'
+      double precision energia(0:n_ene)
+      include 'COMMON.FFIELD'
+      include 'COMMON.DERIV'
+      include 'COMMON.INTERACT'
+      include 'COMMON.SBRIDGE'
+      include 'COMMON.CHAIN'
+      include 'COMMON.VAR'
+      include 'COMMON.LOCAL'
+
+c      write(iout,'(a,i2)')'Calling etotal_short ipot=',ipot
+c      call flush(iout)
+      if (modecalc.eq.12.or.modecalc.eq.14) then
+#ifdef MPI
+        if (fg_rank.eq.0) call int_from_cart1(.false.)
+#else
+        call int_from_cart1(.false.)
+#endif
+      endif
+#ifdef MPI      
+c      write(iout,*) "ETOTAL_SHORT Processor",fg_rank,
+c     & " absolute rank",myrank," nfgtasks",nfgtasks
+c      call flush(iout)
+      if (nfgtasks.gt.1) then
+        time00=MPI_Wtime()
+C FG slaves call the following matching MPI_Bcast in ERGASTULUM
+        if (fg_rank.eq.0) then
+          call MPI_Bcast(2,1,MPI_INTEGER,king,FG_COMM,IERROR)
+c          write (iout,*) "Processor",myrank," BROADCAST iorder"
+c          call flush(iout)
+C FG master sets up the WEIGHTS_ array which will be broadcast to the 
+C FG slaves as WEIGHTS array.
+          weights_(1)=wsc
+          weights_(2)=wscp
+          weights_(3)=welec
+          weights_(4)=wcorr
+          weights_(5)=wcorr5
+          weights_(6)=wcorr6
+          weights_(7)=wel_loc
+          weights_(8)=wturn3
+          weights_(9)=wturn4
+          weights_(10)=wturn6
+          weights_(11)=wang
+          weights_(12)=wscloc
+          weights_(13)=wtor
+          weights_(14)=wtor_d
+          weights_(15)=wstrain
+          weights_(16)=wvdwpp
+          weights_(17)=wbond
+          weights_(18)=scal14
+          weights_(21)=wsccor
+C FG Master broadcasts the WEIGHTS_ array
+          call MPI_Bcast(weights_(1),n_ene,
+     &        MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
+        else
+C FG slaves receive the WEIGHTS array
+          call MPI_Bcast(weights(1),n_ene,
+     &        MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
+          wsc=weights(1)
+          wscp=weights(2)
+          welec=weights(3)
+          wcorr=weights(4)
+          wcorr5=weights(5)
+          wcorr6=weights(6)
+          wel_loc=weights(7)
+          wturn3=weights(8)
+          wturn4=weights(9)
+          wturn6=weights(10)
+          wang=weights(11)
+          wscloc=weights(12)
+          wtor=weights(13)
+          wtor_d=weights(14)
+          wstrain=weights(15)
+          wvdwpp=weights(16)
+          wbond=weights(17)
+          scal14=weights(18)
+          wsccor=weights(21)
+        endif
+c        write (iout,*),"Processor",myrank," BROADCAST weights"
+        call MPI_Bcast(c(1,1),maxres6,MPI_DOUBLE_PRECISION,
+     &    king,FG_COMM,IERR)
+c        write (iout,*) "Processor",myrank," BROADCAST c"
+        call MPI_Bcast(dc(1,1),maxres6,MPI_DOUBLE_PRECISION,
+     &    king,FG_COMM,IERR)
+c        write (iout,*) "Processor",myrank," BROADCAST dc"
+        call MPI_Bcast(dc_norm(1,1),maxres6,MPI_DOUBLE_PRECISION,
+     &    king,FG_COMM,IERR)
+c        write (iout,*) "Processor",myrank," BROADCAST dc_norm"
+        call MPI_Bcast(theta(1),nres,MPI_DOUBLE_PRECISION,
+     &    king,FG_COMM,IERR)
+c        write (iout,*) "Processor",myrank," BROADCAST theta"
+        call MPI_Bcast(phi(1),nres,MPI_DOUBLE_PRECISION,
+     &    king,FG_COMM,IERR)
+c        write (iout,*) "Processor",myrank," BROADCAST phi"
+        call MPI_Bcast(alph(1),nres,MPI_DOUBLE_PRECISION,
+     &    king,FG_COMM,IERR)
+c        write (iout,*) "Processor",myrank," BROADCAST alph"
+        call MPI_Bcast(omeg(1),nres,MPI_DOUBLE_PRECISION,
+     &    king,FG_COMM,IERR)
+c        write (iout,*) "Processor",myrank," BROADCAST omeg"
+        call MPI_Bcast(vbld(1),2*nres,MPI_DOUBLE_PRECISION,
+     &    king,FG_COMM,IERR)
+c        write (iout,*) "Processor",myrank," BROADCAST vbld"
+        call MPI_Bcast(vbld_inv(1),2*nres,MPI_DOUBLE_PRECISION,
+     &    king,FG_COMM,IERR)
+         time_Bcast=time_Bcast+MPI_Wtime()-time00
+c        write (iout,*) "Processor",myrank," BROADCAST vbld_inv"
+      endif
+c      write (iout,*) 'Processor',myrank,
+c     &  ' calling etotal_short ipot=',ipot
+c      call flush(iout)
+c      print *,'Processor',myrank,' nnt=',nnt,' nct=',nct
+#endif     
+c      call int_from_cart1(.false.)
+C
+C Compute the side-chain and electrostatic interaction energy
+C
+      goto (101,102,103,104,105,106) ipot
+C Lennard-Jones potential.
+  101 call elj_short(evdw)
+cd    print '(a)','Exit ELJ'
+      goto 107
+C Lennard-Jones-Kihara potential (shifted).
+  102 call eljk_short(evdw)
+      goto 107
+C Berne-Pechukas potential (dilated LJ, angular dependence).
+  103 call ebp_short(evdw)
+      goto 107
+C Gay-Berne potential (shifted LJ, angular dependence).
+  104 call egb_short(evdw,evdw_p,evdw_m)
+      goto 107
+C Gay-Berne-Vorobjev potential (shifted LJ, angular dependence).
+  105 call egbv_short(evdw)
+      goto 107
+C Soft-sphere potential - already dealt with in the long-range part
+  106 evdw=0.0d0
+c  106 call e_softsphere_short(evdw)
+C
+C Calculate electrostatic (H-bonding) energy of the main chain.
+C
+  107 continue
+c
+c Calculate the short-range part of Evdwpp
+c
+      call evdwpp_short(evdw1)
+c
+c Calculate the short-range part of ESCp
+c
+      if (ipot.lt.6) then
+        call escp_short(evdw2,evdw2_14)
+      endif
+c
+c Calculate the bond-stretching energy
+c
+      call ebond(estr)
+C 
+C Calculate the disulfide-bridge and other energy and the contributions
+C from other distance constraints.
+      call edis(ehpb)
+C
+C Calculate the virtual-bond-angle energy.
+C
+      call ebend(ebe)
+C
+C Calculate the SC local energy.
+C
+      call vec_and_deriv
+      call esc(escloc)
+C
+C Calculate the virtual-bond torsional energy.
+C
+      call etor(etors,edihcnstr)
+C
+C 6/23/01 Calculate double-torsional energy
+C
+      call etor_d(etors_d)
+C
+C 21/5/07 Calculate local sicdechain correlation energy
+C
+      if (wsccor.gt.0.0d0) then
+        call eback_sc_corr(esccor)
+      else
+        esccor=0.0d0
+      endif
+C
+C Put energy components into an array
+C
+      do i=1,n_ene
+        energia(i)=0.0d0
+      enddo
+      energia(1)=evdw
+#ifdef SCP14
+      energia(2)=evdw2-evdw2_14
+      energia(18)=evdw2_14
+#else
+      energia(2)=evdw2
+      energia(18)=0.0d0
+#endif
+#ifdef SPLITELE
+      energia(16)=evdw1
+#else
+      energia(3)=evdw1
+#endif
+      energia(11)=ebe
+      energia(12)=escloc
+      energia(13)=etors
+      energia(14)=etors_d
+      energia(15)=ehpb
+      energia(17)=estr
+      energia(19)=edihcnstr
+      energia(21)=esccor
+      energia(22)=evdw_p
+      energia(23)=evdw_m
+c      write (iout,*) "ETOTAL_SHORT before SUM_ENERGY"
+      call flush(iout)
+      call sum_energy(energia,.true.)
+c      write (iout,*) "Exit ETOTAL_SHORT"
+      call flush(iout)
+      return
+      end