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