source/unres/src-HCD-5D/energy_split-sep.F

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