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