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