Added homology restraints modified from Pawel and Magda's code
[unres.git] / source / unres / src_MD-restraints / energy_split-sep.F
1       subroutine etotal_long(energia)
2       implicit real*8 (a-h,o-z)
3       include 'DIMENSIONS'
4 c
5 c Compute the long-range slow-varying contributions to the energy
6 c
7 #ifndef ISNAN
8       external proc_proc
9 #ifdef WINPGI
10 cMS$ATTRIBUTES C ::  proc_proc
11 #endif
12 #endif
13 #ifdef MPI
14       include "mpif.h"
15       double precision weights_(n_ene)
16 #endif
17       include 'COMMON.SETUP'
18       include 'COMMON.IOUNITS'
19       double precision energia(0:n_ene)
20       include 'COMMON.FFIELD'
21       include 'COMMON.DERIV'
22       include 'COMMON.INTERACT'
23       include 'COMMON.SBRIDGE'
24       include 'COMMON.CHAIN'
25       include 'COMMON.VAR'
26       include 'COMMON.LOCAL'
27       include 'COMMON.MD'
28 c      write(iout,'(a,i2)')'Calling etotal_long ipot=',ipot
29       if (modecalc.eq.12.or.modecalc.eq.14) then
30 #ifdef MPI
31 c        if (fg_rank.eq.0) call int_from_cart1(.false.)
32 #else
33         call int_from_cart1(.false.)
34 #endif
35       endif
36 #ifdef MPI      
37 c      write(iout,*) "ETOTAL_LONG Processor",fg_rank,
38 c     & " absolute rank",myrank," nfgtasks",nfgtasks
39       call flush(iout)
40       if (nfgtasks.gt.1) then
41         time00=MPI_Wtime()
42 C FG slaves call the following matching MPI_Bcast in ERGASTULUM
43         if (fg_rank.eq.0) then
44           call MPI_Bcast(3,1,MPI_INTEGER,king,FG_COMM,IERROR)
45 c          write (iout,*) "Processor",myrank," BROADCAST iorder"
46 c          call flush(iout)
47 C FG master sets up the WEIGHTS_ array which will be broadcast to the 
48 C FG slaves as WEIGHTS array.
49           weights_(1)=wsc
50           weights_(2)=wscp
51           weights_(3)=welec
52           weights_(4)=wcorr
53           weights_(5)=wcorr5
54           weights_(6)=wcorr6
55           weights_(7)=wel_loc
56           weights_(8)=wturn3
57           weights_(9)=wturn4
58           weights_(10)=wturn6
59           weights_(11)=wang
60           weights_(12)=wscloc
61           weights_(13)=wtor
62           weights_(14)=wtor_d
63           weights_(15)=wstrain
64           weights_(16)=wvdwpp
65           weights_(17)=wbond
66           weights_(18)=scal14
67           weights_(21)=wsccor
68 C FG Master broadcasts the WEIGHTS_ array
69           call MPI_Bcast(weights_(1),n_ene,
70      &        MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
71         else
72 C FG slaves receive the WEIGHTS array
73           call MPI_Bcast(weights(1),n_ene,
74      &        MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
75           wsc=weights(1)
76           wscp=weights(2)
77           welec=weights(3)
78           wcorr=weights(4)
79           wcorr5=weights(5)
80           wcorr6=weights(6)
81           wel_loc=weights(7)
82           wturn3=weights(8)
83           wturn4=weights(9)
84           wturn6=weights(10)
85           wang=weights(11)
86           wscloc=weights(12)
87           wtor=weights(13)
88           wtor_d=weights(14)
89           wstrain=weights(15)
90           wvdwpp=weights(16)
91           wbond=weights(17)
92           scal14=weights(18)
93           wsccor=weights(21)
94         endif
95         call MPI_Bcast(dc(1,1),6*nres,MPI_DOUBLE_PRECISION,
96      &    king,FG_COMM,IERR)
97          time_Bcast=time_Bcast+MPI_Wtime()-time00
98          time_Bcastw=time_Bcastw+MPI_Wtime()-time00
99 c        call chainbuild_cart
100 c        call int_from_cart1(.false.)
101       endif
102 c      write (iout,*) 'Processor',myrank,
103 c     &  ' calling etotal_short ipot=',ipot
104 c      call flush(iout)
105 c      print *,'Processor',myrank,' nnt=',nnt,' nct=',nct
106 #endif     
107 cd    print *,'nnt=',nnt,' nct=',nct
108 C
109 C Compute the side-chain and electrostatic interaction energy
110 C
111       goto (101,102,103,104,105,106) ipot
112 C Lennard-Jones potential.
113   101 call elj_long(evdw)
114 cd    print '(a)','Exit ELJ'
115       goto 107
116 C Lennard-Jones-Kihara potential (shifted).
117   102 call eljk_long(evdw)
118       goto 107
119 C Berne-Pechukas potential (dilated LJ, angular dependence).
120   103 call ebp_long(evdw)
121       goto 107
122 C Gay-Berne potential (shifted LJ, angular dependence).
123   104 call egb_long(evdw,evdw_p,evdw_m)
124       goto 107
125 C Gay-Berne-Vorobjev potential (shifted LJ, angular dependence).
126   105 call egbv_long(evdw)
127       goto 107
128 C Soft-sphere potential
129   106 call e_softsphere(evdw)
130 C
131 C Calculate electrostatic (H-bonding) energy of the main chain.
132 C
133   107 continue
134       call vec_and_deriv
135       if (ipot.lt.6) then
136 #ifdef SPLITELE
137          if (welec.gt.0d0.or.wvdwpp.gt.0d0.or.wel_loc.gt.0d0.or.
138      &       wturn3.gt.0d0.or.wturn4.gt.0d0 .or. wcorr.gt.0.0d0
139      &       .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.d0
140      &       .or. wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0 ) then
141 #else
142          if (welec.gt.0d0.or.wel_loc.gt.0d0.or.
143      &       wturn3.gt.0d0.or.wturn4.gt.0d0 .or. wcorr.gt.0.0d0
144      &       .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.d0
145      &       .or. wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0 ) then
146 #endif
147            call eelec_scale(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
148          else
149             ees=0
150             evdw1=0
151             eel_loc=0
152             eello_turn3=0
153             eello_turn4=0
154          endif
155       else
156 c        write (iout,*) "Soft-spheer ELEC potential"
157         call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3,
158      &   eello_turn4)
159       endif
160 C
161 C Calculate excluded-volume interaction energy between peptide groups
162 C and side chains.
163 C
164       if (ipot.lt.6) then
165        if(wscp.gt.0d0) then
166         call escp_long(evdw2,evdw2_14)
167        else
168         evdw2=0
169         evdw2_14=0
170        endif
171       else
172         call escp_soft_sphere(evdw2,evdw2_14)
173       endif
174
175 C 12/1/95 Multi-body terms
176 C
177       n_corr=0
178       n_corr1=0
179       if ((wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0 
180      &    .or. wturn6.gt.0.0d0) .and. ipot.lt.6) then
181          call multibody_eello(ecorr,ecorr5,ecorr6,eturn6,n_corr,n_corr1)
182 c         write (2,*) 'n_corr=',n_corr,' n_corr1=',n_corr1,
183 c     &" ecorr",ecorr," ecorr5",ecorr5," ecorr6",ecorr6," eturn6",eturn6
184       else
185          ecorr=0.0d0
186          ecorr5=0.0d0
187          ecorr6=0.0d0
188          eturn6=0.0d0
189       endif
190       if ((wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) .and. ipot.lt.6) then
191          call multibody_hb(ecorr,ecorr5,ecorr6,n_corr,n_corr1)
192       endif
193
194 C If performing constraint dynamics, call the constraint energy
195 C  after the equilibration time
196       if(usampl.and.totT.gt.eq_time) then
197          call EconstrQ   
198          call Econstr_back
199       else
200          Uconst=0.0d0
201          Uconst_back=0.0d0
202       endif
203
204 C Sum the energies
205 C
206       do i=1,n_ene
207         energia(i)=0.0d0
208       enddo
209       energia(1)=evdw
210 #ifdef SCP14
211       energia(2)=evdw2-evdw2_14
212       energia(18)=evdw2_14
213 #else
214       energia(2)=evdw2
215       energia(18)=0.0d0
216 #endif
217 #ifdef SPLITELE
218       energia(3)=ees
219       energia(16)=evdw1
220 #else
221       energia(3)=ees+evdw1
222       energia(16)=0.0d0
223 #endif
224       energia(4)=ecorr
225       energia(5)=ecorr5
226       energia(6)=ecorr6
227       energia(7)=eel_loc
228       energia(8)=eello_turn3
229       energia(9)=eello_turn4
230       energia(10)=eturn6
231       energia(20)=Uconst+Uconst_back
232       energia(22)=evdw_p
233       energia(23)=evdw_m
234       call sum_energy(energia,.true.)
235 c      write (iout,*) "Exit ETOTAL_LONG"
236       call flush(iout)
237       return
238       end
239 c------------------------------------------------------------------------------
240       subroutine etotal_short(energia)
241       implicit real*8 (a-h,o-z)
242       include 'DIMENSIONS'
243 c
244 c Compute the short-range fast-varying contributions to the energy
245 c
246 #ifndef ISNAN
247       external proc_proc
248 #ifdef WINPGI
249 cMS$ATTRIBUTES C ::  proc_proc
250 #endif
251 #endif
252 #ifdef MPI
253       include "mpif.h"
254       double precision weights_(n_ene)
255 #endif
256       include 'COMMON.SETUP'
257       include 'COMMON.IOUNITS'
258       double precision energia(0:n_ene)
259       include 'COMMON.FFIELD'
260       include 'COMMON.DERIV'
261       include 'COMMON.INTERACT'
262       include 'COMMON.SBRIDGE'
263       include 'COMMON.CHAIN'
264       include 'COMMON.VAR'
265       include 'COMMON.LOCAL'
266       include 'COMMON.CONTROL'
267
268 c      write(iout,'(a,i2)')'Calling etotal_short ipot=',ipot
269 c      call flush(iout)
270       if (modecalc.eq.12.or.modecalc.eq.14) then
271 #ifdef MPI
272         if (fg_rank.eq.0) call int_from_cart1(.false.)
273 #else
274         call int_from_cart1(.false.)
275 #endif
276       endif
277 #ifdef MPI      
278 c      write(iout,*) "ETOTAL_SHORT Processor",fg_rank,
279 c     & " absolute rank",myrank," nfgtasks",nfgtasks
280 c      call flush(iout)
281       if (nfgtasks.gt.1) then
282         time00=MPI_Wtime()
283 C FG slaves call the following matching MPI_Bcast in ERGASTULUM
284         if (fg_rank.eq.0) then
285           call MPI_Bcast(2,1,MPI_INTEGER,king,FG_COMM,IERROR)
286 c          write (iout,*) "Processor",myrank," BROADCAST iorder"
287 c          call flush(iout)
288 C FG master sets up the WEIGHTS_ array which will be broadcast to the 
289 C FG slaves as WEIGHTS array.
290           weights_(1)=wsc
291           weights_(2)=wscp
292           weights_(3)=welec
293           weights_(4)=wcorr
294           weights_(5)=wcorr5
295           weights_(6)=wcorr6
296           weights_(7)=wel_loc
297           weights_(8)=wturn3
298           weights_(9)=wturn4
299           weights_(10)=wturn6
300           weights_(11)=wang
301           weights_(12)=wscloc
302           weights_(13)=wtor
303           weights_(14)=wtor_d
304           weights_(15)=wstrain
305           weights_(16)=wvdwpp
306           weights_(17)=wbond
307           weights_(18)=scal14
308           weights_(21)=wsccor
309 C FG Master broadcasts the WEIGHTS_ array
310           call MPI_Bcast(weights_(1),n_ene,
311      &        MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
312         else
313 C FG slaves receive the WEIGHTS array
314           call MPI_Bcast(weights(1),n_ene,
315      &        MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
316           wsc=weights(1)
317           wscp=weights(2)
318           welec=weights(3)
319           wcorr=weights(4)
320           wcorr5=weights(5)
321           wcorr6=weights(6)
322           wel_loc=weights(7)
323           wturn3=weights(8)
324           wturn4=weights(9)
325           wturn6=weights(10)
326           wang=weights(11)
327           wscloc=weights(12)
328           wtor=weights(13)
329           wtor_d=weights(14)
330           wstrain=weights(15)
331           wvdwpp=weights(16)
332           wbond=weights(17)
333           scal14=weights(18)
334           wsccor=weights(21)
335         endif
336 c        write (iout,*),"Processor",myrank," BROADCAST weights"
337         call MPI_Bcast(c(1,1),maxres6,MPI_DOUBLE_PRECISION,
338      &    king,FG_COMM,IERR)
339 c        write (iout,*) "Processor",myrank," BROADCAST c"
340         call MPI_Bcast(dc(1,1),maxres6,MPI_DOUBLE_PRECISION,
341      &    king,FG_COMM,IERR)
342 c        write (iout,*) "Processor",myrank," BROADCAST dc"
343         call MPI_Bcast(dc_norm(1,1),maxres6,MPI_DOUBLE_PRECISION,
344      &    king,FG_COMM,IERR)
345 c        write (iout,*) "Processor",myrank," BROADCAST dc_norm"
346         call MPI_Bcast(theta(1),nres,MPI_DOUBLE_PRECISION,
347      &    king,FG_COMM,IERR)
348 c        write (iout,*) "Processor",myrank," BROADCAST theta"
349         call MPI_Bcast(phi(1),nres,MPI_DOUBLE_PRECISION,
350      &    king,FG_COMM,IERR)
351 c        write (iout,*) "Processor",myrank," BROADCAST phi"
352         call MPI_Bcast(alph(1),nres,MPI_DOUBLE_PRECISION,
353      &    king,FG_COMM,IERR)
354 c        write (iout,*) "Processor",myrank," BROADCAST alph"
355         call MPI_Bcast(omeg(1),nres,MPI_DOUBLE_PRECISION,
356      &    king,FG_COMM,IERR)
357 c        write (iout,*) "Processor",myrank," BROADCAST omeg"
358         call MPI_Bcast(vbld(1),2*nres,MPI_DOUBLE_PRECISION,
359      &    king,FG_COMM,IERR)
360 c        write (iout,*) "Processor",myrank," BROADCAST vbld"
361         call MPI_Bcast(vbld_inv(1),2*nres,MPI_DOUBLE_PRECISION,
362      &    king,FG_COMM,IERR)
363          time_Bcast=time_Bcast+MPI_Wtime()-time00
364 c        write (iout,*) "Processor",myrank," BROADCAST vbld_inv"
365       endif
366 c      write (iout,*) 'Processor',myrank,
367 c     &  ' calling etotal_short ipot=',ipot
368 c      call flush(iout)
369 c      print *,'Processor',myrank,' nnt=',nnt,' nct=',nct
370 #endif     
371 c      call int_from_cart1(.false.)
372 C
373 C Compute the side-chain and electrostatic interaction energy
374 C
375       goto (101,102,103,104,105,106) ipot
376 C Lennard-Jones potential.
377   101 call elj_short(evdw)
378 cd    print '(a)','Exit ELJ'
379       goto 107
380 C Lennard-Jones-Kihara potential (shifted).
381   102 call eljk_short(evdw)
382       goto 107
383 C Berne-Pechukas potential (dilated LJ, angular dependence).
384   103 call ebp_short(evdw)
385       goto 107
386 C Gay-Berne potential (shifted LJ, angular dependence).
387   104 call egb_short(evdw,evdw_p,evdw_m)
388       goto 107
389 C Gay-Berne-Vorobjev potential (shifted LJ, angular dependence).
390   105 call egbv_short(evdw)
391       goto 107
392 C Soft-sphere potential - already dealt with in the long-range part
393   106 evdw=0.0d0
394 c  106 call e_softsphere_short(evdw)
395 C
396 C Calculate electrostatic (H-bonding) energy of the main chain.
397 C
398   107 continue
399 c
400 c Calculate the short-range part of Evdwpp
401 c
402       call evdwpp_short(evdw1)
403 c
404 c Calculate the short-range part of ESCp
405 c
406       if (ipot.lt.6) then
407         call escp_short(evdw2,evdw2_14)
408       endif
409 c
410 c Calculate the bond-stretching energy
411 c
412       call ebond(estr)
413
414 C Calculate the disulfide-bridge and other energy and the contributions
415 C from other distance constraints.
416       call edis(ehpb)
417 C
418 C Calculate the virtual-bond-angle energy.
419 C
420       call ebend(ebe)
421 C
422 C Calculate the SC local energy.
423 C
424       call vec_and_deriv
425       call esc(escloc)
426 C
427 C Calculate the virtual-bond torsional energy.
428 C
429       call etor(etors,edihcnstr)
430 c
431 c Homology restraints
432 c
433       if (constr_homology.ge.1) then
434         call e_modeller(ehomology_constr)
435       else
436         ehomology_constr=0.0d0
437       endif
438 C
439 C 6/23/01 Calculate double-torsional energy
440 C
441       call etor_d(etors_d)
442 C
443 C 21/5/07 Calculate local sicdechain correlation energy
444 C
445       if (wsccor.gt.0.0d0) then
446         call eback_sc_corr(esccor)
447       else
448         esccor=0.0d0
449       endif
450 C
451 C Put energy components into an array
452 C
453       do i=1,n_ene
454         energia(i)=0.0d0
455       enddo
456       energia(1)=evdw
457 #ifdef SCP14
458       energia(2)=evdw2-evdw2_14
459       energia(18)=evdw2_14
460 #else
461       energia(2)=evdw2
462       energia(18)=0.0d0
463 #endif
464 #ifdef SPLITELE
465       energia(16)=evdw1
466 #else
467       energia(3)=evdw1
468 #endif
469       energia(11)=ebe
470       energia(12)=escloc
471       energia(13)=etors
472       energia(14)=etors_d
473       energia(15)=ehpb
474       energia(17)=estr
475       energia(19)=edihcnstr
476       energia(21)=esccor
477       energia(22)=evdw_p
478       energia(23)=evdw_m
479       energia(24)=ehomology_constr
480 c      write (iout,*) "ETOTAL_SHORT before SUM_ENERGY"
481       call flush(iout)
482       call sum_energy(energia,.true.)
483 c      write (iout,*) "Exit ETOTAL_SHORT"
484       call flush(iout)
485       return
486       end