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