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