source/wham/src-M/differ

   1 1c1
   2 <       subroutine etotal(energia,fact)
   3 ---
   4 >       subroutine etotal(energia)
   5 4,5d3
   6 <       include 'DIMENSIONS.ZSCOPT'
   7 <
   8 8d5
   9 < #endif
  10 12,18d8
  11 <
  12 <       include 'COMMON.IOUNITS'
  13 <       double precision energia(0:max_ene),energia1(0:max_ene+1)
  14 < #ifdef MPL
  15 <       include 'COMMON.INFO'
  16 <       external d_vadd
  17 <       integer ready
  18 19a10,17
  19 > #ifdef MPI
  20 >       include "mpif.h"
  21 >       double precision weights_(n_ene)
  22 > #endif
  23 >       include 'COMMON.SETUP'
  24 >       include 'COMMON.IOUNITS'
  25 >       double precision energia(0:n_ene)
  26 >       include 'COMMON.LOCAL'
  27 25,28c23,98
  28 <       double precision fact(6)
  29 < cd      write(iout, '(a,i2)')'Calling etotal ipot=',ipot
  30 < cd    print *,'nnt=',nnt,' nct=',nct
  31 < C
  32 ---
  33 >       include 'COMMON.VAR'
  34 >       include 'COMMON.MD'
  35 >       include 'COMMON.CONTROL'
  36 >       include 'COMMON.TIME1'
  37 > #ifdef MPI
  38 > c      print*,"ETOTAL Processor",fg_rank," absolute rank",myrank,
  39 > c     & " nfgtasks",nfgtasks
  40 >       if (nfgtasks.gt.1) then
  41 >         time00=MPI_Wtime()
  42 > C FG slaves call the following matching MPI_Bcast in ERGASTULUM
  43 >         if (fg_rank.eq.0) then
  44 >           call MPI_Bcast(0,1,MPI_INTEGER,king,FG_COMM,IERROR)
  45 > c          print *,"Processor",myrank," BROADCAST iorder"
  46 > C FG master sets up the WEIGHTS_ array which will be broadcast to the
  47 > C FG slaves as WEIGHTS array.
  48 >           weights_(1)=wsc
  49 >           weights_(2)=wscp
  50 >           weights_(3)=welec
  51 >           weights_(4)=wcorr
  52 >           weights_(5)=wcorr5
  53 >           weights_(6)=wcorr6
  54 >           weights_(7)=wel_loc
  55 >           weights_(8)=wturn3
  56 >           weights_(9)=wturn4
  57 >           weights_(10)=wturn6
  58 >           weights_(11)=wang
  59 >           weights_(12)=wscloc
  60 >           weights_(13)=wtor
  61 >           weights_(14)=wtor_d
  62 >           weights_(15)=wstrain
  63 >           weights_(16)=wvdwpp
  64 >           weights_(17)=wbond
  65 >           weights_(18)=scal14
  66 >           weights_(21)=wsccor
  67 > C FG Master broadcasts the WEIGHTS_ array
  68 >           call MPI_Bcast(weights_(1),n_ene,
  69 >      &        MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
  70 >         else
  71 > C FG slaves receive the WEIGHTS array
  72 >           call MPI_Bcast(weights(1),n_ene,
  73 >      &        MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
  74 >           wsc=weights(1)
  75 >           wscp=weights(2)
  76 >           welec=weights(3)
  77 >           wcorr=weights(4)
  78 >           wcorr5=weights(5)
  79 >           wcorr6=weights(6)
  80 >           wel_loc=weights(7)
  81 >           wturn3=weights(8)
  82 >           wturn4=weights(9)
  83 >           wturn6=weights(10)
  84 >           wang=weights(11)
  85 >           wscloc=weights(12)
  86 >           wtor=weights(13)
  87 >           wtor_d=weights(14)
  88 >           wstrain=weights(15)
  89 >           wvdwpp=weights(16)
  90 >           wbond=weights(17)
  91 >           scal14=weights(18)
  92 >           wsccor=weights(21)
  93 >         endif
  94 >         time_Bcast=time_Bcast+MPI_Wtime()-time00
  95 >         time_Bcastw=time_Bcastw+MPI_Wtime()-time00
  96 > c        call chainbuild_cart
  97 >       endif
  98 > c      print *,'Processor',myrank,' calling etotal ipot=',ipot
  99 > c      print *,'Processor',myrank,' nnt=',nnt,' nct=',nct
 100 > #else
 101 > c      if (modecalc.eq.12.or.modecalc.eq.14) then
 102 > c        call int_from_cart1(.false.)
 103 > c      endif
 104 > #endif
 105 > #ifdef TIMING
 106 >       time00=MPI_Wtime()
 107 > #endif
 108 > C
 109 31c101
 110 <       goto (101,102,103,104,105) ipot
 111 ---
 112 >       goto (101,102,103,104,105,106) ipot
 113 33c103
 114 <   101 call elj(evdw,evdw_t)
 115 ---
 116 >   101 call elj(evdw)
 117 35c105
 118 <       goto 106
 119 ---
 120 >       goto 107
 121 37,38c107,108
 122 <   102 call eljk(evdw,evdw_t)
 123 <       goto 106
 124 ---
 125 >   102 call eljk(evdw)
 126 >       goto 107
 127 40,41c110,111
 128 <   103 call ebp(evdw,evdw_t)
 129 <       goto 106
 130 ---
 131 >   103 call ebp(evdw)
 132 >       goto 107
 133 43,44c113,114
 134 <   104 call egb(evdw,evdw_t)
 135 <       goto 106
 136 ---
 137 >   104 call egb(evdw)
 138 >       goto 107
 139 46c116,119
 140 <   105 call egbv(evdw,evdw_t)
 141 ---
 142 >   105 call egbv(evdw)
 143 >       goto 107
 144 > C Soft-sphere potential
 145 >   106 call e_softsphere(evdw)
 146 50c123,158
 147 <   106 call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
 148 ---
 149 >   107 continue
 150 > c      print *,"Processor",myrank," computed USCSC"
 151 > #ifdef TIMING
 152 >       time01=MPI_Wtime()
 153 > #endif
 154 >       call vec_and_deriv
 155 > #ifdef TIMING
 156 >       time_vec=time_vec+MPI_Wtime()-time01
 157 > #endif
 158 > c      print *,"Processor",myrank," left VEC_AND_DERIV"
 159 >       if (ipot.lt.6) then
 160 > #ifdef SPLITELE
 161 >          if (welec.gt.0d0.or.wvdwpp.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 > #else
 166 >          if (welec.gt.0d0.or.wel_loc.gt.0d0.or.
 167 >      &       wturn3.gt.0d0.or.wturn4.gt.0d0 .or. wcorr.gt.0.0d0
 168 >      &       .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.d0
 169 >      &       .or. wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0 ) then
 170 > #endif
 171 >             call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
 172 >          else
 173 >             ees=0.0d0
 174 >             evdw1=0.0d0
 175 >             eel_loc=0.0d0
 176 >             eello_turn3=0.0d0
 177 >             eello_turn4=0.0d0
 178 >          endif
 179 >       else
 180 > c        write (iout,*) "Soft-spheer ELEC potential"
 181 >         call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3,
 182 >      &   eello_turn4)
 183 >       endif
 184 > c      print *,"Processor",myrank," computed UELEC"
 185 55c163,173
 186 <       call escp(evdw2,evdw2_14)
 187 ---
 188 >       if (ipot.lt.6) then
 189 >        if(wscp.gt.0d0) then
 190 >         call escp(evdw2,evdw2_14)
 191 >        else
 192 >         evdw2=0
 193 >         evdw2_14=0
 194 >        endif
 195 >       else
 196 > c        write (iout,*) "Soft-sphere SCP potential"
 197 >         call escp_soft_sphere(evdw2,evdw2_14)
 198 >       endif
 199 60d177
 200 < c      write (iout,*) "estr",estr
 201 70,71c187,192
 202 <       call ebend(ebe)
 203 < cd    print *,'Bend energy finished.'
 204 ---
 205 >       if (wang.gt.0d0) then
 206 >         call ebend(ebe)
 207 >       else
 208 >         ebe=0
 209 >       endif
 210 > c      print *,"Processor",myrank," computed UB"
 211 76c197
 212 < cd    print *,'SCLOC energy finished.'
 213 ---
 214 > c      print *,"Processor",myrank," computed USC"
 215 81c202,208
 216 <       call etor(etors,edihcnstr,fact(1))
 217 ---
 218 >       if (wtor.gt.0) then
 219 >        call etor(etors,edihcnstr)
 220 >       else
 221 >        etors=0
 222 >        edihcnstr=0
 223 >       endif
 224 > c      print *,"Processor",myrank," computed Utor"
 225 85c212,217
 226 <       call etor_d(etors_d,fact(2))
 227 ---
 228 >       if (wtor_d.gt.0) then
 229 >        call etor_d(etors_d)
 230 >       else
 231 >        etors_d=0
 232 >       endif
 233 > c      print *,"Processor",myrank," computed Utord"
 234 89c221,226
 235 <       call eback_sc_corr(esccor)
 236 ---
 237 >       if (wsccor.gt.0.0d0) then
 238 >         call eback_sc_corr(esccor)
 239 >       else
 240 >         esccor=0.0d0
 241 >       endif
 242 > c      print *,"Processor",myrank," computed Usccorr"
 243 95,97c232,233
 244 <       if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0
 245 <      &    .or. wturn6.gt.0.0d0) then
 246 < c         print *,"calling multibody_eello"
 247 ---
 248 >       if ((wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0
 249 >      &    .or. wturn6.gt.0.0d0) .and. ipot.lt.6) then
 250 99,100c235,241
 251 < c         write (*,*) 'n_corr=',n_corr,' n_corr1=',n_corr1
 252 < c         print *,ecorr,ecorr5,ecorr6,eturn6
 253 ---
 254 > cd         write(2,*)'multibody_eello n_corr=',n_corr,' n_corr1=',n_corr1,
 255 > cd     &" ecorr",ecorr," ecorr5",ecorr5," ecorr6",ecorr6," eturn6",eturn6
 256 >       else
 257 >          ecorr=0.0d0
 258 >          ecorr5=0.0d0
 259 >          ecorr6=0.0d0
 260 >          eturn6=0.0d0
 261 102c243
 262 <       if (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) then
 263 ---
 264 >       if ((wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) .and. ipot.lt.6) then
 265 103a245
 266 > cd         write (iout,*) "multibody_hb ecorr",ecorr
 267 105,123c247,259
 268 < c      write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t
 269 < #ifdef SPLITELE
 270 <       etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2+welec*fact(1)*ees
 271 <      & +wvdwpp*evdw1
 272 <      & +wang*ebe+wtor*fact(1)*etors+wscloc*escloc
 273 <      & +wstrain*ehpb+nss*ebr+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
 274 <      & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4
 275 <      & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
 276 <      & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
 277 <      & +wbond*estr+wsccor*fact(1)*esccor
 278 < #else
 279 <       etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2
 280 <      & +welec*fact(1)*(ees+evdw1)
 281 <      & +wang*ebe+wtor*fact(1)*etors+wscloc*escloc
 282 <      & +wstrain*ehpb+nss*ebr+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
 283 <      & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4
 284 <      & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
 285 <      & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
 286 <      & +wbond*estr+wsccor*fact(1)*esccor
 287 ---
 288 > c      print *,"Processor",myrank," computed Ucorr"
 289 > C
 290 > C If performing constraint dynamics, call the constraint energy
 291 > C  after the equilibration time
 292 >       if(usampl.and.totT.gt.eq_time) then
 293 >          call EconstrQ
 294 >          call Econstr_back
 295 >       else
 296 >          Uconst=0.0d0
 297 >          Uconst_back=0.0d0
 298 >       endif
 299 > #ifdef TIMING
 300 >       time_enecalc=time_enecalc+MPI_Wtime()-time00
 301 125c261,267
 302 <       energia(0)=etot
 303 ---
 304 > c      print *,"Processor",myrank," computed Uconstr"
 305 > #ifdef TIMING
 306 >       time00=MPI_Wtime()
 307 > #endif
 308 > c
 309 > C Sum the energies
 310 > C
 311 129c271
 312 <       energia(17)=evdw2_14
 313 ---
 314 >       energia(18)=evdw2_14
 315 132c274
 316 <       energia(17)=0.0d0
 317 ---
 318 >       energia(18)=0.0d0
 319 153,156c295,402
 320 <       energia(18)=estr
 321 <       energia(19)=esccor
 322 <       energia(20)=edihcnstr
 323 <       energia(21)=evdw_t
 324 ---
 325 >       energia(19)=edihcnstr
 326 >       energia(17)=estr
 327 >       energia(20)=Uconst+Uconst_back
 328 >       energia(21)=esccor
 329 > c      print *," Processor",myrank," calls SUM_ENERGY"
 330 >       call sum_energy(energia,.true.)
 331 > c      print *," Processor",myrank," left SUM_ENERGY"
 332 > #ifdef TIMING
 333 >       time_sumene=time_sumene+MPI_Wtime()-time00
 334 > #endif
 335 >       return
 336 >       end
 337 > c-------------------------------------------------------------------------------
 338 >       subroutine sum_energy(energia,reduce)
 339 >       implicit real*8 (a-h,o-z)
 340 >       include 'DIMENSIONS'
 341 > #ifndef ISNAN
 342 >       external proc_proc
 343 > #ifdef WINPGI
 344 > cMS$ATTRIBUTES C ::  proc_proc
 345 > #endif
 346 > #endif
 347 > #ifdef MPI
 348 >       include "mpif.h"
 349 > #endif
 350 >       include 'COMMON.SETUP'
 351 >       include 'COMMON.IOUNITS'
 352 >       double precision energia(0:n_ene),enebuff(0:n_ene+1)
 353 >       include 'COMMON.FFIELD'
 354 >       include 'COMMON.DERIV'
 355 >       include 'COMMON.INTERACT'
 356 >       include 'COMMON.SBRIDGE'
 357 >       include 'COMMON.CHAIN'
 358 >       include 'COMMON.VAR'
 359 >       include 'COMMON.CONTROL'
 360 >       include 'COMMON.TIME1'
 361 >       logical reduce
 362 > #ifdef MPI
 363 >       if (nfgtasks.gt.1 .and. reduce) then
 364 > #ifdef DEBUG
 365 >         write (iout,*) "energies before REDUCE"
 366 >         call enerprint(energia)
 367 >         call flush(iout)
 368 > #endif
 369 >         do i=0,n_ene
 370 >           enebuff(i)=energia(i)
 371 >         enddo
 372 >         time00=MPI_Wtime()
 373 >         call MPI_Barrier(FG_COMM,IERR)
 374 >         time_barrier_e=time_barrier_e+MPI_Wtime()-time00
 375 >         time00=MPI_Wtime()
 376 >         call MPI_Reduce(enebuff(0),energia(0),n_ene+1,
 377 >      &    MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
 378 > #ifdef DEBUG
 379 >         write (iout,*) "energies after REDUCE"
 380 >         call enerprint(energia)
 381 >         call flush(iout)
 382 > #endif
 383 >         time_Reduce=time_Reduce+MPI_Wtime()-time00
 384 >       endif
 385 >       if (fg_rank.eq.0) then
 386 > #endif
 387 >       evdw=energia(1)
 388 > #ifdef SCP14
 389 >       evdw2=energia(2)+energia(18)
 390 >       evdw2_14=energia(18)
 391 > #else
 392 >       evdw2=energia(2)
 393 > #endif
 394 > #ifdef SPLITELE
 395 >       ees=energia(3)
 396 >       evdw1=energia(16)
 397 > #else
 398 >       ees=energia(3)
 399 >       evdw1=0.0d0
 400 > #endif
 401 >       ecorr=energia(4)
 402 >       ecorr5=energia(5)
 403 >       ecorr6=energia(6)
 404 >       eel_loc=energia(7)
 405 >       eello_turn3=energia(8)
 406 >       eello_turn4=energia(9)
 407 >       eturn6=energia(10)
 408 >       ebe=energia(11)
 409 >       escloc=energia(12)
 410 >       etors=energia(13)
 411 >       etors_d=energia(14)
 412 >       ehpb=energia(15)
 413 >       edihcnstr=energia(19)
 414 >       estr=energia(17)
 415 >       Uconst=energia(20)
 416 >       esccor=energia(21)
 417 > #ifdef SPLITELE
 418 >       etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1
 419 >      & +wang*ebe+wtor*etors+wscloc*escloc
 420 >      & +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5
 421 >      & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
 422 >      & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
 423 >      & +wbond*estr+Uconst+wsccor*esccor
 424 > #else
 425 >       etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1)
 426 >      & +wang*ebe+wtor*etors+wscloc*escloc
 427 >      & +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5
 428 >      & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
 429 >      & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
 430 >      & +wbond*estr+Uconst+wsccor*esccor
 431 > #endif
 432 >       energia(0)=etot
 433 173,174c419,464
 434 < #ifdef MPL
 435 < c     endif
 436 ---
 437 > #ifdef MPI
 438 >       endif
 439 > #endif
 440 >       return
 441 >       end
 442 > c-------------------------------------------------------------------------------
 443 >       subroutine sum_gradient
 444 >       implicit real*8 (a-h,o-z)
 445 >       include 'DIMENSIONS'
 446 > #ifndef ISNAN
 447 >       external proc_proc
 448 > #ifdef WINPGI
 449 > cMS$ATTRIBUTES C ::  proc_proc
 450 > #endif
 451 > #endif
 452 > #ifdef MPI
 453 >       include 'mpif.h'
 454 >       double precision gradbufc(3,maxres),gradbufx(3,maxres),
 455 >      &  glocbuf(4*maxres),gradbufc_sum(3,maxres)
 456 > #endif
 457 >       include 'COMMON.SETUP'
 458 >       include 'COMMON.IOUNITS'
 459 >       include 'COMMON.FFIELD'
 460 >       include 'COMMON.DERIV'
 461 >       include 'COMMON.INTERACT'
 462 >       include 'COMMON.SBRIDGE'
 463 >       include 'COMMON.CHAIN'
 464 >       include 'COMMON.VAR'
 465 >       include 'COMMON.CONTROL'
 466 >       include 'COMMON.TIME1'
 467 >       include 'COMMON.MAXGRAD'
 468 > #ifdef TIMING
 469 >       time01=MPI_Wtime()
 470 > #endif
 471 > #ifdef DEBUG
 472 >       write (iout,*) "sum_gradient gvdwc, gvdwx"
 473 >       do i=1,nres
 474 >         write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)')
 475 >      &   i,(gvdwx(j,i),j=1,3),(gvdwc(j,i),j=1,3)
 476 >       enddo
 477 >       call flush(iout)
 478 > #endif
 479 > #ifdef MPI
 480 > C FG slaves call the following matching MPI_Bcast in ERGASTULUM
 481 >         if (nfgtasks.gt.1 .and. fg_rank.eq.0)
 482 >      &    call MPI_Bcast(1,1,MPI_INTEGER,king,FG_COMM,IERROR)
 483 176d465
 484 <       if (calc_grad) then
 485 178c467,468
 486 < C Sum up the components of the Cartesian gradient.
 487 ---
 488 > C 9/29/08 AL Transform parts of gradients in site coordinates to the gradient
 489 > C            in virtual-bond-vector coordinates
 490 179a470,488
 491 > #ifdef DEBUG
 492 > c      write (iout,*) "gel_loc gel_loc_long and gel_loc_loc"
 493 > c      do i=1,nres-1
 494 > c        write (iout,'(i5,3f10.5,2x,3f10.5,2x,f10.5)')
 495 > c     &   i,(gel_loc(j,i),j=1,3),(gel_loc_long(j,i),j=1,3),gel_loc_loc(i)
 496 > c      enddo
 497 > c      write (iout,*) "gel_loc_tur3 gel_loc_turn4"
 498 > c      do i=1,nres-1
 499 > c        write (iout,'(i5,3f10.5,2x,f10.5)')
 500 > c     &  i,(gcorr4_turn(j,i),j=1,3),gel_loc_turn4(i)
 501 > c      enddo
 502 >       write (iout,*) "gradcorr5 gradcorr5_long gradcorr5_loc"
 503 >       do i=1,nres
 504 >         write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)')
 505 >      &   i,(gradcorr5(j,i),j=1,3),(gradcorr5_long(j,i),j=1,3),
 506 >      &   g_corr5_loc(i)
 507 >       enddo
 508 >       call flush(iout)
 509 > #endif
 510 183,198c492,500
 511 <           gradc(j,i,icg)=wsc*gvdwc(j,i)+wscp*gvdwc_scp(j,i)+
 512 <      &                welec*fact(1)*gelc(j,i)+wvdwpp*gvdwpp(j,i)+
 513 <      &                wbond*gradb(j,i)+
 514 <      &                wstrain*ghpbc(j,i)+
 515 <      &                wcorr*fact(3)*gradcorr(j,i)+
 516 <      &                wel_loc*fact(2)*gel_loc(j,i)+
 517 <      &                wturn3*fact(2)*gcorr3_turn(j,i)+
 518 <      &                wturn4*fact(3)*gcorr4_turn(j,i)+
 519 <      &                wcorr5*fact(4)*gradcorr5(j,i)+
 520 <      &                wcorr6*fact(5)*gradcorr6(j,i)+
 521 <      &                wturn6*fact(5)*gcorr6_turn(j,i)+
 522 <      &                wsccor*fact(2)*gsccorc(j,i)
 523 <           gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+
 524 <      &                  wbond*gradbx(j,i)+
 525 <      &                  wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+
 526 <      &                  wsccor*fact(2)*gsccorx(j,i)
 527 ---
 528 >           gradbufc(j,i)=wsc*gvdwc(j,i)+
 529 >      &                wscp*(gvdwc_scp(j,i)+gvdwc_scpp(j,i))+
 530 >      &                welec*gelc_long(j,i)+wvdwpp*gvdwpp(j,i)+
 531 >      &                wel_loc*gel_loc_long(j,i)+
 532 >      &                wcorr*gradcorr_long(j,i)+
 533 >      &                wcorr5*gradcorr5_long(j,i)+
 534 >      &                wcorr6*gradcorr6_long(j,i)+
 535 >      &                wturn6*gcorr6_turn_long(j,i)+
 536 >      &                wstrain*ghpbc(j,i)
 537 199a502
 538 >       enddo
 539 203,204c506,508
 540 <           gradc(j,i,icg)=wsc*gvdwc(j,i)+wscp*gvdwc_scp(j,i)+
 541 <      &                welec*fact(1)*gelc(j,i)+wstrain*ghpbc(j,i)+
 542 ---
 543 >           gradbufc(j,i)=wsc*gvdwc(j,i)+
 544 >      &                wscp*(gvdwc_scp(j,i)+gvdwc_scpp(j,i))+
 545 >      &                welec*gelc_long(j,i)+
 546 206,213c510,670
 547 <      &                wcorr*fact(3)*gradcorr(j,i)+
 548 <      &                wel_loc*fact(2)*gel_loc(j,i)+
 549 <      &                wturn3*fact(2)*gcorr3_turn(j,i)+
 550 <      &                wturn4*fact(3)*gcorr4_turn(j,i)+
 551 <      &                wcorr5*fact(4)*gradcorr5(j,i)+
 552 <      &                wcorr6*fact(5)*gradcorr6(j,i)+
 553 <      &                wturn6*fact(5)*gcorr6_turn(j,i)+
 554 <      &                wsccor*fact(2)*gsccorc(j,i)
 555 ---
 556 >      &                wel_loc*gel_loc_long(j,i)+
 557 >      &                wcorr*gradcorr_long(j,i)+
 558 >      &                wcorr5*gradcorr5_long(j,i)+
 559 >      &                wcorr6*gradcorr6_long(j,i)+
 560 >      &                wturn6*gcorr6_turn_long(j,i)+
 561 >      &                wstrain*ghpbc(j,i)
 562 >         enddo
 563 >       enddo
 564 > #endif
 565 > #ifdef MPI
 566 >       if (nfgtasks.gt.1) then
 567 >       time00=MPI_Wtime()
 568 > #ifdef DEBUG
 569 >       write (iout,*) "gradbufc before allreduce"
 570 >       do i=1,nres
 571 >         write (iout,'(i3,3f10.5)') i,(gradbufc(j,i),j=1,3)
 572 >       enddo
 573 >       call flush(iout)
 574 > #endif
 575 >       do i=1,nres
 576 >         do j=1,3
 577 >           gradbufc_sum(j,i)=gradbufc(j,i)
 578 >         enddo
 579 >       enddo
 580 > c      call MPI_AllReduce(gradbufc(1,1),gradbufc_sum(1,1),3*nres,
 581 > c     &    MPI_DOUBLE_PRECISION,MPI_SUM,FG_COMM,IERR)
 582 > c      time_reduce=time_reduce+MPI_Wtime()-time00
 583 > #ifdef DEBUG
 584 > c      write (iout,*) "gradbufc_sum after allreduce"
 585 > c      do i=1,nres
 586 > c        write (iout,'(i3,3f10.5)') i,(gradbufc_sum(j,i),j=1,3)
 587 > c      enddo
 588 > c      call flush(iout)
 589 > #endif
 590 > #ifdef TIMING
 591 > c      time_allreduce=time_allreduce+MPI_Wtime()-time00
 592 > #endif
 593 >       do i=nnt,nres
 594 >         do k=1,3
 595 >           gradbufc(k,i)=0.0d0
 596 >         enddo
 597 >       enddo
 598 > #ifdef DEBUG
 599 >       write (iout,*) "igrad_start",igrad_start," igrad_end",igrad_end
 600 >       write (iout,*) (i," jgrad_start",jgrad_start(i),
 601 >      &                  " jgrad_end  ",jgrad_end(i),
 602 >      &                  i=igrad_start,igrad_end)
 603 > #endif
 604 > c
 605 > c Obsolete and inefficient code; we can make the effort O(n) and, therefore,
 606 > c do not parallelize this part.
 607 > c
 608 > c      do i=igrad_start,igrad_end
 609 > c        do j=jgrad_start(i),jgrad_end(i)
 610 > c          do k=1,3
 611 > c            gradbufc(k,i)=gradbufc(k,i)+gradbufc_sum(k,j)
 612 > c          enddo
 613 > c        enddo
 614 > c      enddo
 615 >       do j=1,3
 616 >         gradbufc(j,nres-1)=gradbufc_sum(j,nres)
 617 >       enddo
 618 >       do i=nres-2,nnt,-1
 619 >         do j=1,3
 620 >           gradbufc(j,i)=gradbufc(j,i+1)+gradbufc_sum(j,i+1)
 621 >         enddo
 622 >       enddo
 623 > #ifdef DEBUG
 624 >       write (iout,*) "gradbufc after summing"
 625 >       do i=1,nres
 626 >         write (iout,'(i3,3f10.5)') i,(gradbufc(j,i),j=1,3)
 627 >       enddo
 628 >       call flush(iout)
 629 > #endif
 630 >       else
 631 > #endif
 632 > #ifdef DEBUG
 633 >       write (iout,*) "gradbufc"
 634 >       do i=1,nres
 635 >         write (iout,'(i3,3f10.5)') i,(gradbufc(j,i),j=1,3)
 636 >       enddo
 637 >       call flush(iout)
 638 > #endif
 639 >       do i=1,nres
 640 >         do j=1,3
 641 >           gradbufc_sum(j,i)=gradbufc(j,i)
 642 >           gradbufc(j,i)=0.0d0
 643 >         enddo
 644 >       enddo
 645 >       do j=1,3
 646 >         gradbufc(j,nres-1)=gradbufc_sum(j,nres)
 647 >       enddo
 648 >       do i=nres-2,nnt,-1
 649 >         do j=1,3
 650 >           gradbufc(j,i)=gradbufc(j,i+1)+gradbufc_sum(j,i+1)
 651 >         enddo
 652 >       enddo
 653 > c      do i=nnt,nres-1
 654 > c        do k=1,3
 655 > c          gradbufc(k,i)=0.0d0
 656 > c        enddo
 657 > c        do j=i+1,nres
 658 > c          do k=1,3
 659 > c            gradbufc(k,i)=gradbufc(k,i)+gradbufc(k,j)
 660 > c          enddo
 661 > c        enddo
 662 > c      enddo
 663 > #ifdef DEBUG
 664 >       write (iout,*) "gradbufc after summing"
 665 >       do i=1,nres
 666 >         write (iout,'(i3,3f10.5)') i,(gradbufc(j,i),j=1,3)
 667 >       enddo
 668 >       call flush(iout)
 669 > #endif
 670 > #ifdef MPI
 671 >       endif
 672 > #endif
 673 >       do k=1,3
 674 >         gradbufc(k,nres)=0.0d0
 675 >       enddo
 676 >       do i=1,nct
 677 >         do j=1,3
 678 > #ifdef SPLITELE
 679 >           gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+
 680 >      &                wel_loc*gel_loc(j,i)+
 681 >      &                0.5d0*(wscp*gvdwc_scpp(j,i)+
 682 >      &                welec*gelc_long(j,i)+wvdwpp*gvdwpp(j,i)+
 683 >      &                wel_loc*gel_loc_long(j,i)+
 684 >      &                wcorr*gradcorr_long(j,i)+
 685 >      &                wcorr5*gradcorr5_long(j,i)+
 686 >      &                wcorr6*gradcorr6_long(j,i)+
 687 >      &                wturn6*gcorr6_turn_long(j,i))+
 688 >      &                wbond*gradb(j,i)+
 689 >      &                wcorr*gradcorr(j,i)+
 690 >      &                wturn3*gcorr3_turn(j,i)+
 691 >      &                wturn4*gcorr4_turn(j,i)+
 692 >      &                wcorr5*gradcorr5(j,i)+
 693 >      &                wcorr6*gradcorr6(j,i)+
 694 >      &                wturn6*gcorr6_turn(j,i)+
 695 >      &                wsccor*gsccorc(j,i)
 696 >      &               +wscloc*gscloc(j,i)
 697 > #else
 698 >           gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+
 699 >      &                wel_loc*gel_loc(j,i)+
 700 >      &                0.5d0*(wscp*gvdwc_scpp(j,i)+
 701 >      &                welec*gelc_long(j,i)
 702 >      &                wel_loc*gel_loc_long(j,i)+
 703 >      &                wcorr*gcorr_long(j,i)+
 704 >      &                wcorr5*gradcorr5_long(j,i)+
 705 >      &                wcorr6*gradcorr6_long(j,i)+
 706 >      &                wturn6*gcorr6_turn_long(j,i))+
 707 >      &                wbond*gradb(j,i)+
 708 >      &                wcorr*gradcorr(j,i)+
 709 >      &                wturn3*gcorr3_turn(j,i)+
 710 >      &                wturn4*gcorr4_turn(j,i)+
 711 >      &                wcorr5*gradcorr5(j,i)+
 712 >      &                wcorr6*gradcorr6(j,i)+
 713 >      &                wturn6*gcorr6_turn(j,i)+
 714 >      &                wsccor*gsccorc(j,i)
 715 >      &               +wscloc*gscloc(j,i)
 716 > #endif
 717 217c674,675
 718 <      &                  wsccor*fact(1)*gsccorx(j,i)
 719 ---
 720 >      &                  wsccor*gsccorx(j,i)
 721 >      &                 +wscloc*gsclocx(j,i)
 722 219c677,681
 723 < #endif
 724 ---
 725 >       enddo
 726 > #ifdef DEBUG
 727 >       write (iout,*) "gloc before adding corr"
 728 >       do i=1,4*nres
 729 >         write (iout,*) i,gloc(i,icg)
 730 221,222c683
 731 <
 732 <
 733 ---
 734 > #endif
 735 224,231c685,697
 736 <         gloc(i,icg)=gloc(i,icg)+wcorr*fact(3)*gcorr_loc(i)
 737 <      &   +wcorr5*fact(4)*g_corr5_loc(i)
 738 <      &   +wcorr6*fact(5)*g_corr6_loc(i)
 739 <      &   +wturn4*fact(3)*gel_loc_turn4(i)
 740 <      &   +wturn3*fact(2)*gel_loc_turn3(i)
 741 <      &   +wturn6*fact(5)*gel_loc_turn6(i)
 742 <      &   +wel_loc*fact(2)*gel_loc_loc(i)
 743 <      &   +wsccor*fact(1)*gsccor_loc(i)
 744 ---
 745 >         gloc(i,icg)=gloc(i,icg)+wcorr*gcorr_loc(i)
 746 >      &   +wcorr5*g_corr5_loc(i)
 747 >      &   +wcorr6*g_corr6_loc(i)
 748 >      &   +wturn4*gel_loc_turn4(i)
 749 >      &   +wturn3*gel_loc_turn3(i)
 750 >      &   +wturn6*gel_loc_turn6(i)
 751 >      &   +wel_loc*gel_loc_loc(i)
 752 >      &   +wsccor*gsccor_loc(i)
 753 >       enddo
 754 > #ifdef DEBUG
 755 >       write (iout,*) "gloc after adding corr"
 756 >       do i=1,4*nres
 757 >         write (iout,*) i,gloc(i,icg)
 758 232a699,727
 759 > #endif
 760 > #ifdef MPI
 761 >       if (nfgtasks.gt.1) then
 762 >         do j=1,3
 763 >           do i=1,nres
 764 >             gradbufc(j,i)=gradc(j,i,icg)
 765 >             gradbufx(j,i)=gradx(j,i,icg)
 766 >           enddo
 767 >         enddo
 768 >         do i=1,4*nres
 769 >           glocbuf(i)=gloc(i,icg)
 770 >         enddo
 771 >         time00=MPI_Wtime()
 772 >         call MPI_Barrier(FG_COMM,IERR)
 773 >         time_barrier_g=time_barrier_g+MPI_Wtime()-time00
 774 >         time00=MPI_Wtime()
 775 >         call MPI_Reduce(gradbufc(1,1),gradc(1,1,icg),3*nres,
 776 >      &    MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
 777 >         call MPI_Reduce(gradbufx(1,1),gradx(1,1,icg),3*nres,
 778 >      &    MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
 779 >         call MPI_Reduce(glocbuf(1),gloc(1,icg),4*nres,
 780 >      &    MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
 781 >         time_reduce=time_reduce+MPI_Wtime()-time00
 782 > #ifdef DEBUG
 783 >       write (iout,*) "gloc after reduce"
 784 >       do i=1,4*nres
 785 >         write (iout,*) i,gloc(i,icg)
 786 >       enddo
 787 > #endif
 788 233a729,897
 789 > #endif
 790 >       if (gnorm_check) then
 791 > c
 792 > c Compute the maximum elements of the gradient
 793 > c
 794 >       gvdwc_max=0.0d0
 795 >       gvdwc_scp_max=0.0d0
 796 >       gelc_max=0.0d0
 797 >       gvdwpp_max=0.0d0
 798 >       gradb_max=0.0d0
 799 >       ghpbc_max=0.0d0
 800 >       gradcorr_max=0.0d0
 801 >       gel_loc_max=0.0d0
 802 >       gcorr3_turn_max=0.0d0
 803 >       gcorr4_turn_max=0.0d0
 804 >       gradcorr5_max=0.0d0
 805 >       gradcorr6_max=0.0d0
 806 >       gcorr6_turn_max=0.0d0
 807 >       gsccorc_max=0.0d0
 808 >       gscloc_max=0.0d0
 809 >       gvdwx_max=0.0d0
 810 >       gradx_scp_max=0.0d0
 811 >       ghpbx_max=0.0d0
 812 >       gradxorr_max=0.0d0
 813 >       gsccorx_max=0.0d0
 814 >       gsclocx_max=0.0d0
 815 >       do i=1,nct
 816 >         gvdwc_norm=dsqrt(scalar(gvdwc(1,i),gvdwc(1,i)))
 817 >         if (gvdwc_norm.gt.gvdwc_max) gvdwc_max=gvdwc_norm
 818 >         gvdwc_scp_norm=dsqrt(scalar(gvdwc_scp(1,i),gvdwc_scp(1,i)))
 819 >         if (gvdwc_scp_norm.gt.gvdwc_scp_max)
 820 >      &   gvdwc_scp_max=gvdwc_scp_norm
 821 >         gelc_norm=dsqrt(scalar(gelc(1,i),gelc(1,i)))
 822 >         if (gelc_norm.gt.gelc_max) gelc_max=gelc_norm
 823 >         gvdwpp_norm=dsqrt(scalar(gvdwpp(1,i),gvdwpp(1,i)))
 824 >         if (gvdwpp_norm.gt.gvdwpp_max) gvdwpp_max=gvdwpp_norm
 825 >         gradb_norm=dsqrt(scalar(gradb(1,i),gradb(1,i)))
 826 >         if (gradb_norm.gt.gradb_max) gradb_max=gradb_norm
 827 >         ghpbc_norm=dsqrt(scalar(ghpbc(1,i),ghpbc(1,i)))
 828 >         if (ghpbc_norm.gt.ghpbc_max) ghpbc_max=ghpbc_norm
 829 >         gradcorr_norm=dsqrt(scalar(gradcorr(1,i),gradcorr(1,i)))
 830 >         if (gradcorr_norm.gt.gradcorr_max) gradcorr_max=gradcorr_norm
 831 >         gel_loc_norm=dsqrt(scalar(gel_loc(1,i),gel_loc(1,i)))
 832 >         if (gel_loc_norm.gt.gel_loc_max) gel_loc_max=gel_loc_norm
 833 >         gcorr3_turn_norm=dsqrt(scalar(gcorr3_turn(1,i),
 834 >      &    gcorr3_turn(1,i)))
 835 >         if (gcorr3_turn_norm.gt.gcorr3_turn_max)
 836 >      &    gcorr3_turn_max=gcorr3_turn_norm
 837 >         gcorr4_turn_norm=dsqrt(scalar(gcorr4_turn(1,i),
 838 >      &    gcorr4_turn(1,i)))
 839 >         if (gcorr4_turn_norm.gt.gcorr4_turn_max)
 840 >      &    gcorr4_turn_max=gcorr4_turn_norm
 841 >         gradcorr5_norm=dsqrt(scalar(gradcorr5(1,i),gradcorr5(1,i)))
 842 >         if (gradcorr5_norm.gt.gradcorr5_max)
 843 >      &    gradcorr5_max=gradcorr5_norm
 844 >         gradcorr6_norm=dsqrt(scalar(gradcorr6(1,i),gradcorr6(1,i)))
 845 >         if (gradcorr6_norm.gt.gradcorr6_max) gcorr6_max=gradcorr6_norm
 846 >         gcorr6_turn_norm=dsqrt(scalar(gcorr6_turn(1,i),
 847 >      &    gcorr6_turn(1,i)))
 848 >         if (gcorr6_turn_norm.gt.gcorr6_turn_max)
 849 >      &    gcorr6_turn_max=gcorr6_turn_norm
 850 >         gsccorr_norm=dsqrt(scalar(gsccorc(1,i),gsccorc(1,i)))
 851 >         if (gsccorr_norm.gt.gsccorr_max) gsccorr_max=gsccorr_norm
 852 >         gscloc_norm=dsqrt(scalar(gscloc(1,i),gscloc(1,i)))
 853 >         if (gscloc_norm.gt.gscloc_max) gscloc_max=gscloc_norm
 854 >         gvdwx_norm=dsqrt(scalar(gvdwx(1,i),gvdwx(1,i)))
 855 >         if (gvdwx_norm.gt.gvdwx_max) gvdwx_max=gvdwx_norm
 856 >         gradx_scp_norm=dsqrt(scalar(gradx_scp(1,i),gradx_scp(1,i)))
 857 >         if (gradx_scp_norm.gt.gradx_scp_max)
 858 >      &    gradx_scp_max=gradx_scp_norm
 859 >         ghpbx_norm=dsqrt(scalar(ghpbx(1,i),ghpbx(1,i)))
 860 >         if (ghpbx_norm.gt.ghpbx_max) ghpbx_max=ghpbx_norm
 861 >         gradxorr_norm=dsqrt(scalar(gradxorr(1,i),gradxorr(1,i)))
 862 >         if (gradxorr_norm.gt.gradxorr_max) gradxorr_max=gradxorr_norm
 863 >         gsccorrx_norm=dsqrt(scalar(gsccorx(1,i),gsccorx(1,i)))
 864 >         if (gsccorrx_norm.gt.gsccorrx_max) gsccorrx_max=gsccorrx_norm
 865 >         gsclocx_norm=dsqrt(scalar(gsclocx(1,i),gsclocx(1,i)))
 866 >         if (gsclocx_norm.gt.gsclocx_max) gsclocx_max=gsclocx_norm
 867 >       enddo
 868 >       if (gradout) then
 869 > #ifdef AIX
 870 >         open(istat,file=statname,position="append")
 871 > #else
 872 >         open(istat,file=statname,access="append")
 873 > #endif
 874 >         write (istat,'(1h#,21f10.2)') gvdwc_max,gvdwc_scp_max,
 875 >      &     gelc_max,gvdwpp_max,gradb_max,ghpbc_max,
 876 >      &     gradcorr_max,gel_loc_max,gcorr3_turn_max,gcorr4_turn_max,
 877 >      &     gradcorr5_max,gradcorr6_max,gcorr6_turn_max,gsccorc_max,
 878 >      &     gscloc_max,gvdwx_max,gradx_scp_max,ghpbx_max,gradxorr_max,
 879 >      &     gsccorx_max,gsclocx_max
 880 >         close(istat)
 881 >         if (gvdwc_max.gt.1.0d4) then
 882 >           write (iout,*) "gvdwc gvdwx gradb gradbx"
 883 >           do i=nnt,nct
 884 >             write(iout,'(i5,4(3f10.2,5x))') i,(gvdwc(j,i),gvdwx(j,i),
 885 >      &        gradb(j,i),gradbx(j,i),j=1,3)
 886 >           enddo
 887 >           call pdbout(0.0d0,'cipiszcze',iout)
 888 >           call flush(iout)
 889 >         endif
 890 >       endif
 891 >       endif
 892 > #ifdef DEBUG
 893 >       write (iout,*) "gradc gradx gloc"
 894 >       do i=1,nres
 895 >         write (iout,'(i5,3f10.5,5x,3f10.5,5x,f10.5)')
 896 >      &   i,(gradc(j,i,icg),j=1,3),(gradx(j,i,icg),j=1,3),gloc(i,icg)
 897 >       enddo
 898 > #endif
 899 > #ifdef TIMING
 900 >       time_sumgradient=time_sumgradient+MPI_Wtime()-time01
 901 > #endif
 902 >       return
 903 >       end
 904 > c-------------------------------------------------------------------------------
 905 >       subroutine rescale_weights(t_bath)
 906 >       implicit real*8 (a-h,o-z)
 907 >       include 'DIMENSIONS'
 908 >       include 'COMMON.IOUNITS'
 909 >       include 'COMMON.FFIELD'
 910 >       include 'COMMON.SBRIDGE'
 911 >       double precision kfac /2.4d0/
 912 >       double precision x,x2,x3,x4,x5,licznik /1.12692801104297249644/
 913 > c      facT=temp0/t_bath
 914 > c      facT=2*temp0/(t_bath+temp0)
 915 >       if (rescale_mode.eq.0) then
 916 >         facT=1.0d0
 917 >         facT2=1.0d0
 918 >         facT3=1.0d0
 919 >         facT4=1.0d0
 920 >         facT5=1.0d0
 921 >       else if (rescale_mode.eq.1) then
 922 >         facT=kfac/(kfac-1.0d0+t_bath/temp0)
 923 >         facT2=kfac**2/(kfac**2-1.0d0+(t_bath/temp0)**2)
 924 >         facT3=kfac**3/(kfac**3-1.0d0+(t_bath/temp0)**3)
 925 >         facT4=kfac**4/(kfac**4-1.0d0+(t_bath/temp0)**4)
 926 >         facT5=kfac**5/(kfac**5-1.0d0+(t_bath/temp0)**5)
 927 >       else if (rescale_mode.eq.2) then
 928 >         x=t_bath/temp0
 929 >         x2=x*x
 930 >         x3=x2*x
 931 >         x4=x3*x
 932 >         x5=x4*x
 933 >         facT=licznik/dlog(dexp(x)+dexp(-x))
 934 >         facT2=licznik/dlog(dexp(x2)+dexp(-x2))
 935 >         facT3=licznik/dlog(dexp(x3)+dexp(-x3))
 936 >         facT4=licznik/dlog(dexp(x4)+dexp(-x4))
 937 >         facT5=licznik/dlog(dexp(x5)+dexp(-x5))
 938 >       else
 939 >         write (iout,*) "Wrong RESCALE_MODE",rescale_mode
 940 >         write (*,*) "Wrong RESCALE_MODE",rescale_mode
 941 > #ifdef MPI
 942 >        call MPI_Finalize(MPI_COMM_WORLD,IERROR)
 943 > #endif
 944 >        stop 555
 945 >       endif
 946 >       welec=weights(3)*fact
 947 >       wcorr=weights(4)*fact3
 948 >       wcorr5=weights(5)*fact4
 949 >       wcorr6=weights(6)*fact5
 950 >       wel_loc=weights(7)*fact2
 951 >       wturn3=weights(8)*fact2
 952 >       wturn4=weights(9)*fact3
 953 >       wturn6=weights(10)*fact5
 954 >       wtor=weights(13)*fact
 955 >       wtor_d=weights(14)*fact2
 956 >       wsccor=weights(21)*fact
 957 >
 958 237c901
 959 <       subroutine enerprint(energia,fact)
 960 ---
 961 >       subroutine enerprint(energia)
 962 240d903
 963 <       include 'DIMENSIONS.ZSCOPT'
 964 244c907,908
 965 <       double precision energia(0:max_ene),fact(6)
 966 ---
 967 >       include 'COMMON.MD'
 968 >       double precision energia(0:n_ene)
 969 246c910,911
 970 <       evdw=energia(1)+fact(6)*energia(21)
 971 ---
 972 >       evdw=energia(1)
 973 >       evdw2=energia(2)
 974 248c913
 975 <       evdw2=energia(2)+energia(17)
 976 ---
 977 >       evdw2=energia(2)+energia(18)
 978 268,270c933,936
 979 <       esccor=energia(19)
 980 <       edihcnstr=energia(20)
 981 <       estr=energia(18)
 982 ---
 983 >       edihcnstr=energia(19)
 984 >       estr=energia(17)
 985 >       Uconst=energia(20)
 986 >       esccor=energia(21)
 987 272,279c938,945
 988 <       write (iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),evdw1,
 989 <      &  wvdwpp,
 990 <      &  estr,wbond,ebe,wang,escloc,wscloc,etors,wtor*fact(1),
 991 <      &  etors_d,wtor_d*fact(2),ehpb,wstrain,
 992 <      &  ecorr,wcorr*fact(3),ecorr5,wcorr5*fact(4),ecorr6,wcorr6*fact(5),
 993 <      &  eel_loc,wel_loc*fact(2),eello_turn3,wturn3*fact(2),
 994 <      &  eello_turn4,wturn4*fact(3),eello_turn6,wturn6*fact(5),
 995 <      &  esccor,wsccor*fact(1),edihcnstr,ebr*nss,etot
 996 ---
 997 >       write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,
 998 >      &  estr,wbond,ebe,wang,
 999 >      &  escloc,wscloc,etors,wtor,etors_d,wtor_d,ehpb,wstrain,
1000 >      &  ecorr,wcorr,
1001 >      &  ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
1002 >      &  eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,
1003 >      &  edihcnstr,ebr*nss,
1004 >      &  Uconst,etot
1005 283c949
1006 <      & 'EES=   ',1pE16.6,' WEIGHT=',1pD16.6,' (p-p elec)'/
1007 ---
1008 >      & 'EES=   ',1pE16.6,' WEIGHT=',1pD16.6,' (p-p)'/
1009 301c967,968
1010 <      & 'ESS=   ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
1011 ---
1012 >      & 'ESS=   ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
1013 >      & 'UCONST= ',1pE16.6,' (Constraint energy)'/
1014 304,310c971,977
1015 <       write (iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),estr,wbond,
1016 <      &  ebe,wang,escloc,wscloc,etors,wtor*fact(1),etors_d,wtor_d*fact2,
1017 <      &  ehpb,wstrain,ecorr,wcorr*fact(3),ecorr5,wcorr5*fact(4),
1018 <      &  ecorr6,wcorr6*fact(5),eel_loc,wel_loc*fact(2),
1019 <      &  eello_turn3,wturn3*fact(2),eello_turn4,wturn4*fact(3),
1020 <      &  eello_turn6,wturn6*fact(5),esccor*fact(1),wsccor,
1021 <      &  edihcnstr,ebr*nss,etot
1022 ---
1023 >       write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,
1024 >      &  estr,wbond,ebe,wang,
1025 >      &  escloc,wscloc,etors,wtor,etors_d,wtor_d,ehpb,wstrain,
1026 >      &  ecorr,wcorr,
1027 >      &  ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
1028 >      &  eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr,
1029 >      &  ebr*nss,Uconst,etot
1030 331c998,999
1031 <      & 'ESS=   ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
1032 ---
1033 >      & 'ESS=   ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
1034 >      & 'UCONST=',1pE16.6,' (Constraint energy)'/
1035 337c1005
1036 <       subroutine elj(evdw,evdw_t)
1037 ---
1038 >       subroutine elj(evdw)
1039 344,345d1011
1040 <       include 'DIMENSIONS.ZSCOPT'
1041 <       include "DIMENSIONS.COMPAR"
1042 354d1019
1043 <       include 'COMMON.ENEPS'
1044 360,367c1025
1045 <       integer icant
1046 <       external icant
1047 < cd    print *,'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
1048 <       do i=1,210
1049 <         do j=1,2
1050 <           eneps_temp(j,i)=0.0d0
1051 <         enddo
1052 <       enddo
1053 ---
1054 > c      write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
1055 369d1026
1056 <       evdw_t=0.0d0
1057 400,402d1056
1058 <             ij=icant(itypi,itypj)
1059 <             eneps_temp(1,ij)=eneps_temp(1,ij)+e1/dabs(eps0ij)
1060 <             eneps_temp(2,ij)=eneps_temp(2,ij)+e2/eps0ij
1061 409,414c1063
1062 <             if (bb(itypi,itypj).gt.0.0d0) then
1063 <               evdw=evdw+evdwij
1064 <             else
1065 <               evdw_t=evdw_t+evdwij
1066 <             endif
1067 <             if (calc_grad) then
1068 ---
1069 >             evdw=evdw+evdwij
1070 424a1074,1075
1071 >               gvdwc(k,i)=gvdwc(k,i)-gg(k)
1072 >               gvdwc(k,j)=gvdwc(k,j)+gg(k)
1073 426,431c1077,1081
1074 <             do k=i,j-1
1075 <               do l=1,3
1076 <                 gvdwc(l,k)=gvdwc(l,k)+gg(l)
1077 <               enddo
1078 <             enddo
1079 <             endif
1080 ---
1081 > cgrad            do k=i,j-1
1082 > cgrad              do l=1,3
1083 > cgrad                gvdwc(l,k)=gvdwc(l,k)+gg(l)
1084 > cgrad              enddo
1085 > cgrad            enddo
1086 493d1142
1087 <       if (calc_grad) then
1088 500d1148
1089 <       endif
1090 513c1161
1091 <       subroutine eljk(evdw,evdw_t)
1092 ---
1093 >       subroutine eljk(evdw)
1094 520,521d1167
1095 <       include 'DIMENSIONS.ZSCOPT'
1096 <       include "DIMENSIONS.COMPAR"
1097 528d1173
1098 <       include 'COMMON.ENEPS'
1099 533,534d1177
1100 <       integer icant
1101 <       external icant
1102 536,540d1178
1103 <       do i=1,210
1104 <         do j=1,2
1105 <           eneps_temp(j,i)=0.0d0
1106 <         enddo
1107 <       enddo
1108 542d1179
1109 <       evdw_t=0.0d0
1110 570,573d1206
1111 <             ij=icant(itypi,itypj)
1112 <             eneps_temp(1,ij)=eneps_temp(1,ij)+(e1+a_augm)
1113 <      &        /dabs(eps(itypi,itypj))
1114 <             eneps_temp(2,ij)=eneps_temp(2,ij)+e2/eps(itypi,itypj)
1115 581,586c1214
1116 <             if (bb(itypi,itypj).gt.0.0d0) then
1117 <               evdw=evdw+evdwij
1118 <             else
1119 <               evdw_t=evdw_t+evdwij
1120 <             endif
1121 <             if (calc_grad) then
1122 ---
1123 >             evdw=evdw+evdwij
1124 596a1225,1226
1125 >               gvdwc(k,i)=gvdwc(k,i)-gg(k)
1126 >               gvdwc(k,j)=gvdwc(k,j)+gg(k)
1127 598,603c1228,1232
1128 <             do k=i,j-1
1129 <               do l=1,3
1130 <                 gvdwc(l,k)=gvdwc(l,k)+gg(l)
1131 <               enddo
1132 <             enddo
1133 <             endif
1134 ---
1135 > cgrad            do k=i,j-1
1136 > cgrad              do l=1,3
1137 > cgrad                gvdwc(l,k)=gvdwc(l,k)+gg(l)
1138 > cgrad              enddo
1139 > cgrad            enddo
1140 607d1235
1141 <       if (calc_grad) then
1142 614d1241
1143 <       endif
1144 618c1245
1145 <       subroutine ebp(evdw,evdw_t)
1146 ---
1147 >       subroutine ebp(evdw)
1148 625,626d1251
1149 <       include 'DIMENSIONS.ZSCOPT'
1150 <       include "DIMENSIONS.COMPAR"
1151 634d1258
1152 <       include 'COMMON.ENEPS'
1153 640,646d1263
1154 <       integer icant
1155 <       external icant
1156 <       do i=1,210
1157 <         do j=1,2
1158 <           eneps_temp(j,i)=0.0d0
1159 <         enddo
1160 <       enddo
1161 648d1264
1162 <       evdw_t=0.0d0
1163 649a1266
1164 >       evdw=0.0D0
1165 665a1283
1166 > c        dsci_inv=dsc_inv(itypi)
1167 674a1293
1168 > c            dscj_inv=dsc_inv(itypj)
1169 719,729c1338
1170 <             ij=icant(itypi,itypj)
1171 <             aux=eps1*eps2rt**2*eps3rt**2
1172 <             eneps_temp(1,ij)=eneps_temp(1,ij)+e1*aux
1173 <      &        /dabs(eps(itypi,itypj))
1174 <             eneps_temp(2,ij)=eneps_temp(2,ij)+e2*aux/eps(itypi,itypj)
1175 <             if (bb(itypi,itypj).gt.0.0d0) then
1176 <               evdw=evdw+evdwij
1177 <             else
1178 <               evdw_t=evdw_t+evdwij
1179 <             endif
1180 <             if (calc_grad) then
1181 ---
1182 >             evdw=evdw+evdwij
1183 752d1360
1184 <             endif
1185 760c1368
1186 <       subroutine egb(evdw,evdw_t)
1187 ---
1188 >       subroutine egb(evdw)
1189 767,768d1374
1190 <       include 'DIMENSIONS.ZSCOPT'
1191 <       include "DIMENSIONS.COMPAR"
1192 776d1381
1193 <       include 'COMMON.ENEPS'
1194 778a1384
1195 >       include 'COMMON.CONTROL'
1196 780,787c1386,1387
1197 <       common /srutu/icall
1198 <       integer icant
1199 <       external icant
1200 <       do i=1,210
1201 <         do j=1,2
1202 <           eneps_temp(j,i)=0.0d0
1203 <         enddo
1204 <       enddo
1205 ---
1206 >       evdw=0.0D0
1207 > ccccc      energy_dec=.false.
1208 790d1389
1209 <       evdw_t=0.0d0
1210 792c1391
1211 < c      if (icall.gt.0) lprn=.true.
1212 ---
1213 > c     if (icall.eq.0) lprn=.false.
1214 803a1403
1215 > c        dsci_inv=dsc_inv(itypi)
1216 804a1405,1406
1217 > c        write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1218 > c        write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1219 812a1415
1220 > c            dscj_inv=dsc_inv(itypj)
1221 813a1417,1419
1222 > c            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1223 > c     &       1.0d0/vbld(j+nres)
1224 > c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1225 840c1446,1448
1226 < c            write (iout,*) i,j,xj,yj,zj
1227 ---
1228 > c            write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1229 > c            write (iout,*) "j",j," dc_norm",
1230 > c     &       dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1231 848a1457,1458
1232 > c for diagnostics; uncomment
1233 > c            rij_shift=1.2*sig0ij
1234 851a1462,1464
1235 > cd              write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1236 > cd     &        restyp(itypi),i,restyp(itypj),j,
1237 > cd     &        rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1238 862a1476,1477
1239 > c            write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1240 > c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1241 864,876c1479
1242 <             if (bb(itypi,itypj).gt.0) then
1243 <               evdw=evdw+evdwij
1244 <             else
1245 <               evdw_t=evdw_t+evdwij
1246 <             endif
1247 <             ij=icant(itypi,itypj)
1248 <             aux=eps1*eps2rt**2*eps3rt**2
1249 <             eneps_temp(1,ij)=eneps_temp(1,ij)+aux*e1
1250 <      &        /dabs(eps(itypi,itypj))
1251 <             eneps_temp(2,ij)=eneps_temp(2,ij)+aux*e2/eps(itypi,itypj)
1252 < c            write (iout,*) "i",i," j",j," itypi",itypi," itypj",itypj,
1253 < c     &         " ij",ij," eneps",aux*e1/dabs(eps(itypi,itypj)),
1254 < c     &         aux*e2/eps(itypi,itypj)
1255 ---
1256 >             evdw=evdw+evdwij
1257 887c1490,1493
1258 <             if (calc_grad) then
1259 ---
1260 >
1261 >             if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
1262 >      &                        'evdw',i,j,evdwij
1263 >
1264 892a1499
1265 > c            fac=0.0d0
1266 899d1505
1267 <             endif
1268 902a1509,1510
1269 > c      write (iout,*) "Number of loop steps in EGB:",ind
1270 > cccc      energy_dec=.false.
1271 906c1514
1272 <       subroutine egbv(evdw,evdw_t)
1273 ---
1274 >       subroutine egbv(evdw)
1275 913,914d1520
1276 <       include 'DIMENSIONS.ZSCOPT'
1277 <       include "DIMENSIONS.COMPAR"
1278 922d1527
1279 <       include 'COMMON.ENEPS'
1280 927,933d1531
1281 <       integer icant
1282 <       external icant
1283 <       do i=1,210
1284 <         do j=1,2
1285 <           eneps_temp(j,i)=0.0d0
1286 <         enddo
1287 <       enddo
1288 935d1532
1289 <       evdw_t=0.0d0
1290 939c1536
1291 < c      if (icall.gt.0) lprn=.true.
1292 ---
1293 > c     if (icall.eq.0) lprn=.true.
1294 950a1548
1295 > c        dsci_inv=dsc_inv(itypi)
1296 959a1558
1297 > c            dscj_inv=dsc_inv(itypj)
1298 1013,1016c1612,1622
1299 <             if (bb(itypi,itypj).gt.0.0d0) then
1300 <               evdw=evdw+evdwij+e_augm
1301 <             else
1302 <               evdw_t=evdw_t+evdwij+e_augm
1303 ---
1304 >             evdw=evdw+evdwij+e_augm
1305 >             if (lprn) then
1306 >             sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1307 >             epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1308 >             write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1309 >      &        restyp(itypi),i,restyp(itypj),j,
1310 >      &        epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1311 >      &        chi1,chi2,chip1,chip2,
1312 >      &        eps1,eps2rt**2,eps3rt**2,
1313 >      &        om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1314 >      &        evdwij+e_augm
1315 1018,1036d1623
1316 <             ij=icant(itypi,itypj)
1317 <             aux=eps1*eps2rt**2*eps3rt**2
1318 <             eneps_temp(1,ij)=eneps_temp(1,ij)+aux*(e1+e_augm)
1319 <      &        /dabs(eps(itypi,itypj))
1320 <             eneps_temp(2,ij)=eneps_temp(2,ij)+aux*e2/eps(itypi,itypj)
1321 < c            eneps_temp(ij)=eneps_temp(ij)
1322 < c     &         +(evdwij+e_augm)/eps(itypi,itypj)
1323 < c            if (lprn) then
1324 < c            sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1325 < c            epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1326 < c            write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1327 < c     &        restyp(itypi),i,restyp(itypj),j,
1328 < c     &        epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1329 < c     &        chi1,chi2,chip1,chip2,
1330 < c     &        eps1,eps2rt**2,eps3rt**2,
1331 < c     &        om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1332 < c     &        evdwij+e_augm
1333 < c            endif
1334 <             if (calc_grad) then
1335 1048d1634
1336 <             endif
1337 1052d1637
1338 <       return
1339 1059a1645
1340 >       include 'COMMON.IOUNITS'
1341 1072a1659,1663
1342 > c diagnostics only
1343 > c      faceps1_inv=om12
1344 > c      eps1=om12
1345 > c      eps1_om12=1.0d0
1346 > c      write (iout,*) "om12",om12," eps1",eps1
1347 1082a1674,1681
1348 > c diagnostics only
1349 > c      sigsq=1.0d0
1350 > c      sigsq_om1=0.0d0
1351 > c      sigsq_om2=0.0d0
1352 > c      sigsq_om12=0.0d0
1353 > c      write (iout,*) "chiom1",chiom1," chiom2",chiom2," chiom12",chiom12
1354 > c      write (iout,*) "faceps1",faceps1," faceps1_inv",faceps1_inv,
1355 > c     &    " eps1",eps1
1356 1089a1689,1690
1357 > c      write (iout,*) "chipom1",chipom1," chipom2",chipom2,
1358 > c     &  " chipom12",chipom12," facp",facp," facp_inv",facp_inv
1359 1098a1700,1702
1360 > c      write (iout,*) "eps2rt",eps2rt," eps3rt",eps3rt
1361 > c      write (iout,*) "eps2rt_om1",eps2rt_om1," eps2rt_om2",eps2rt_om2,
1362 > c     &  " eps2rt_om12",eps2rt_om12
1363 1107d1710
1364 <       include 'DIMENSIONS.ZSCOPT'
1365 1110a1714
1366 >       include 'COMMON.IOUNITS'
1367 1115a1720,1728
1368 > c diagnostics only
1369 > c      eom1=0.0d0
1370 > c      eom2=0.0d0
1371 > c      eom12=evdwij*eps1_om12
1372 > c end diagnostics
1373 > c      write (iout,*) "eps2der",eps2der," eps3der",eps3der,
1374 > c     &  " sigder",sigder
1375 > c      write (iout,*) "eps1_om12",eps1_om12," eps2rt_om12",eps2rt_om12
1376 > c      write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
1377 1122a1736
1378 > c      write (iout,*) "gg",(gg(k),k=1,3)
1379 1129a1744,1747
1380 > c        write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1381 > c     &            +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
1382 > c        write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1383 > c     &            +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
1384 1134,1137c1752,1759
1385 <       do k=i,j-1
1386 <         do l=1,3
1387 <           gvdwc(l,k)=gvdwc(l,k)+gg(l)
1388 <         enddo
1389 ---
1390 > cgrad      do k=i,j-1
1391 > cgrad        do l=1,3
1392 > cgrad          gvdwc(l,k)=gvdwc(l,k)+gg(l)
1393 > cgrad        enddo
1394 > cgrad      enddo
1395 >       do l=1,3
1396 >         gvdwc(l,i)=gvdwc(l,i)-gg(l)
1397 >         gvdwc(l,j)=gvdwc(l,j)+gg(l)
1398 1141,1142c1763,1768
1399 < c------------------------------------------------------------------------------
1400 <       subroutine vec_and_deriv
1401 ---
1402 > C-----------------------------------------------------------------------
1403 >       subroutine e_softsphere(evdw)
1404 > C
1405 > C This subroutine calculates the interaction energy of nonbonded side chains
1406 > C assuming the LJ potential of interaction.
1407 > C
1408 1145,1146c1771
1409 <       include 'DIMENSIONS.ZSCOPT'
1410 <       include 'COMMON.IOUNITS'
1411 ---
1412 >       parameter (accur=1.0d-10)
1413 1151d1775
1414 <       include 'COMMON.VECTORS'
1415 1154,1247c1778,1815
1416 <       dimension uyder(3,3,2),uzder(3,3,2),vbld_inv_temp(2)
1417 < C Compute the local reference systems. For reference system (i), the
1418 < C X-axis points from CA(i) to CA(i+1), the Y axis is in the
1419 < C CA(i)-CA(i+1)-CA(i+2) plane, and the Z axis is perpendicular to this plane.
1420 <       do i=1,nres-1
1421 < c          if (i.eq.nres-1 .or. itel(i+1).eq.0) then
1422 <           if (i.eq.nres-1) then
1423 < C Case of the last full residue
1424 < C Compute the Z-axis
1425 <             call vecpr(dc_norm(1,i),dc_norm(1,i-1),uz(1,i))
1426 <             costh=dcos(pi-theta(nres))
1427 <             fac=1.0d0/dsqrt(1.0d0-costh*costh)
1428 <             do k=1,3
1429 <               uz(k,i)=fac*uz(k,i)
1430 <             enddo
1431 <             if (calc_grad) then
1432 < C Compute the derivatives of uz
1433 <             uzder(1,1,1)= 0.0d0
1434 <             uzder(2,1,1)=-dc_norm(3,i-1)
1435 <             uzder(3,1,1)= dc_norm(2,i-1)
1436 <             uzder(1,2,1)= dc_norm(3,i-1)
1437 <             uzder(2,2,1)= 0.0d0
1438 <             uzder(3,2,1)=-dc_norm(1,i-1)
1439 <             uzder(1,3,1)=-dc_norm(2,i-1)
1440 <             uzder(2,3,1)= dc_norm(1,i-1)
1441 <             uzder(3,3,1)= 0.0d0
1442 <             uzder(1,1,2)= 0.0d0
1443 <             uzder(2,1,2)= dc_norm(3,i)
1444 <             uzder(3,1,2)=-dc_norm(2,i)
1445 <             uzder(1,2,2)=-dc_norm(3,i)
1446 <             uzder(2,2,2)= 0.0d0
1447 <             uzder(3,2,2)= dc_norm(1,i)
1448 <             uzder(1,3,2)= dc_norm(2,i)
1449 <             uzder(2,3,2)=-dc_norm(1,i)
1450 <             uzder(3,3,2)= 0.0d0
1451 <             endif
1452 < C Compute the Y-axis
1453 <             facy=fac
1454 <             do k=1,3
1455 <               uy(k,i)=fac*(dc_norm(k,i-1)-costh*dc_norm(k,i))
1456 <             enddo
1457 <             if (calc_grad) then
1458 < C Compute the derivatives of uy
1459 <             do j=1,3
1460 <               do k=1,3
1461 <                 uyder(k,j,1)=2*dc_norm(k,i-1)*dc_norm(j,i)
1462 <      &                        -dc_norm(k,i)*dc_norm(j,i-1)
1463 <                 uyder(k,j,2)=-dc_norm(j,i)*dc_norm(k,i)
1464 <               enddo
1465 <               uyder(j,j,1)=uyder(j,j,1)-costh
1466 <               uyder(j,j,2)=1.0d0+uyder(j,j,2)
1467 <             enddo
1468 <             do j=1,2
1469 <               do k=1,3
1470 <                 do l=1,3
1471 <                   uygrad(l,k,j,i)=uyder(l,k,j)
1472 <                   uzgrad(l,k,j,i)=uzder(l,k,j)
1473 <                 enddo
1474 <               enddo
1475 <             enddo
1476 <             call unormderiv(uy(1,i),uyder(1,1,1),facy,uygrad(1,1,1,i))
1477 <             call unormderiv(uy(1,i),uyder(1,1,2),facy,uygrad(1,1,2,i))
1478 <             call unormderiv(uz(1,i),uzder(1,1,1),fac,uzgrad(1,1,1,i))
1479 <             call unormderiv(uz(1,i),uzder(1,1,2),fac,uzgrad(1,1,2,i))
1480 <             endif
1481 <           else
1482 < C Other residues
1483 < C Compute the Z-axis
1484 <             call vecpr(dc_norm(1,i),dc_norm(1,i+1),uz(1,i))
1485 <             costh=dcos(pi-theta(i+2))
1486 <             fac=1.0d0/dsqrt(1.0d0-costh*costh)
1487 <             do k=1,3
1488 <               uz(k,i)=fac*uz(k,i)
1489 <             enddo
1490 <             if (calc_grad) then
1491 < C Compute the derivatives of uz
1492 <             uzder(1,1,1)= 0.0d0
1493 <             uzder(2,1,1)=-dc_norm(3,i+1)
1494 <             uzder(3,1,1)= dc_norm(2,i+1)
1495 <             uzder(1,2,1)= dc_norm(3,i+1)
1496 <             uzder(2,2,1)= 0.0d0
1497 <             uzder(3,2,1)=-dc_norm(1,i+1)
1498 <             uzder(1,3,1)=-dc_norm(2,i+1)
1499 <             uzder(2,3,1)= dc_norm(1,i+1)
1500 <             uzder(3,3,1)= 0.0d0
1501 <             uzder(1,1,2)= 0.0d0
1502 <             uzder(2,1,2)= dc_norm(3,i)
1503 <             uzder(3,1,2)=-dc_norm(2,i)
1504 <             uzder(1,2,2)=-dc_norm(3,i)
1505 <             uzder(2,2,2)= 0.0d0
1506 <             uzder(3,2,2)= dc_norm(1,i)
1507 <             uzder(1,3,2)= dc_norm(2,i)
1508 <             uzder(2,3,2)=-dc_norm(1,i)
1509 <             uzder(3,3,2)= 0.0d0
1510 ---
1511 >       include 'COMMON.TORSION'
1512 >       include 'COMMON.SBRIDGE'
1513 >       include 'COMMON.NAMES'
1514 >       include 'COMMON.IOUNITS'
1515 >       include 'COMMON.CONTACTS'
1516 >       dimension gg(3)
1517 > cd    print *,'Entering Esoft_sphere nnt=',nnt,' nct=',nct
1518 >       evdw=0.0D0
1519 >       do i=iatsc_s,iatsc_e
1520 >         itypi=itype(i)
1521 >         if (itypi.eq.21) cycle
1522 >         itypi1=itype(i+1)
1523 >         xi=c(1,nres+i)
1524 >         yi=c(2,nres+i)
1525 >         zi=c(3,nres+i)
1526 > C
1527 > C Calculate SC interaction energy.
1528 > C
1529 >         do iint=1,nint_gr(i)
1530 > cd        write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
1531 > cd   &                  'iend=',iend(i,iint)
1532 >           do j=istart(i,iint),iend(i,iint)
1533 >             itypj=itype(j)
1534 >             if (itypj.eq.21) cycle
1535 >             xj=c(1,nres+j)-xi
1536 >             yj=c(2,nres+j)-yi
1537 >             zj=c(3,nres+j)-zi
1538 >             rij=xj*xj+yj*yj+zj*zj
1539 > c           write (iout,*)'i=',i,' j=',j,' itypi=',itypi,' itypj=',itypj
1540 >             r0ij=r0(itypi,itypj)
1541 >             r0ijsq=r0ij*r0ij
1542 > c            print *,i,j,r0ij,dsqrt(rij)
1543 >             if (rij.lt.r0ijsq) then
1544 >               evdwij=0.25d0*(rij-r0ijsq)**2
1545 >               fac=rij-r0ijsq
1546 >             else
1547 >               evdwij=0.0d0
1548 >               fac=0.0d0
1549 1249,1250c1817,1823
1550 < C Compute the Y-axis
1551 <             facy=fac
1552 ---
1553 >             evdw=evdw+evdwij
1554 > C
1555 > C Calculate the components of the gradient in DC and X
1556 > C
1557 >             gg(1)=xj*fac
1558 >             gg(2)=yj*fac
1559 >             gg(3)=zj*fac
1560 1252,1263c1825,1828
1561 <               uy(k,i)=facy*(dc_norm(k,i+1)-costh*dc_norm(k,i))
1562 <             enddo
1563 <             if (calc_grad) then
1564 < C Compute the derivatives of uy
1565 <             do j=1,3
1566 <               do k=1,3
1567 <                 uyder(k,j,1)=2*dc_norm(k,i+1)*dc_norm(j,i)
1568 <      &                        -dc_norm(k,i)*dc_norm(j,i+1)
1569 <                 uyder(k,j,2)=-dc_norm(j,i)*dc_norm(k,i)
1570 <               enddo
1571 <               uyder(j,j,1)=uyder(j,j,1)-costh
1572 <               uyder(j,j,2)=1.0d0+uyder(j,j,2)
1573 ---
1574 >               gvdwx(k,i)=gvdwx(k,i)-gg(k)
1575 >               gvdwx(k,j)=gvdwx(k,j)+gg(k)
1576 >               gvdwc(k,i)=gvdwc(k,i)-gg(k)
1577 >               gvdwc(k,j)=gvdwc(k,j)+gg(k)
1578 1265,1277c1830,1898
1579 <             do j=1,2
1580 <               do k=1,3
1581 <                 do l=1,3
1582 <                   uygrad(l,k,j,i)=uyder(l,k,j)
1583 <                   uzgrad(l,k,j,i)=uzder(l,k,j)
1584 <                 enddo
1585 <               enddo
1586 <             enddo
1587 <             call unormderiv(uy(1,i),uyder(1,1,1),facy,uygrad(1,1,1,i))
1588 <             call unormderiv(uy(1,i),uyder(1,1,2),facy,uygrad(1,1,2,i))
1589 <             call unormderiv(uz(1,i),uzder(1,1,1),fac,uzgrad(1,1,1,i))
1590 <             call unormderiv(uz(1,i),uzder(1,1,2),fac,uzgrad(1,1,2,i))
1591 <           endif
1592 ---
1593 > cgrad            do k=i,j-1
1594 > cgrad              do l=1,3
1595 > cgrad                gvdwc(l,k)=gvdwc(l,k)+gg(l)
1596 > cgrad              enddo
1597 > cgrad            enddo
1598 >           enddo ! j
1599 >         enddo ! iint
1600 >       enddo ! i
1601 >       return
1602 >       end
1603 > C--------------------------------------------------------------------------
1604 >       subroutine eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3,
1605 >      &              eello_turn4)
1606 > C
1607 > C Soft-sphere potential of p-p interaction
1608 > C
1609 >       implicit real*8 (a-h,o-z)
1610 >       include 'DIMENSIONS'
1611 >       include 'COMMON.CONTROL'
1612 >       include 'COMMON.IOUNITS'
1613 >       include 'COMMON.GEO'
1614 >       include 'COMMON.VAR'
1615 >       include 'COMMON.LOCAL'
1616 >       include 'COMMON.CHAIN'
1617 >       include 'COMMON.DERIV'
1618 >       include 'COMMON.INTERACT'
1619 >       include 'COMMON.CONTACTS'
1620 >       include 'COMMON.TORSION'
1621 >       include 'COMMON.VECTORS'
1622 >       include 'COMMON.FFIELD'
1623 >       dimension ggg(3)
1624 > cd      write(iout,*) 'In EELEC_soft_sphere'
1625 >       ees=0.0D0
1626 >       evdw1=0.0D0
1627 >       eel_loc=0.0d0
1628 >       eello_turn3=0.0d0
1629 >       eello_turn4=0.0d0
1630 >       ind=0
1631 >       do i=iatel_s,iatel_e
1632 >         if (itype(i).eq.21 .or. itype(i+1).eq.21) cycle
1633 >         dxi=dc(1,i)
1634 >         dyi=dc(2,i)
1635 >         dzi=dc(3,i)
1636 >         xmedi=c(1,i)+0.5d0*dxi
1637 >         ymedi=c(2,i)+0.5d0*dyi
1638 >         zmedi=c(3,i)+0.5d0*dzi
1639 >         num_conti=0
1640 > c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
1641 >         do j=ielstart(i),ielend(i)
1642 >           if (itype(j).eq.21 .or. itype(j+1).eq.21) cycle
1643 >           ind=ind+1
1644 >           iteli=itel(i)
1645 >           itelj=itel(j)
1646 >           if (j.eq.i+2 .and. itelj.eq.2) iteli=2
1647 >           r0ij=rpp(iteli,itelj)
1648 >           r0ijsq=r0ij*r0ij
1649 >           dxj=dc(1,j)
1650 >           dyj=dc(2,j)
1651 >           dzj=dc(3,j)
1652 >           xj=c(1,j)+0.5D0*dxj-xmedi
1653 >           yj=c(2,j)+0.5D0*dyj-ymedi
1654 >           zj=c(3,j)+0.5D0*dzj-zmedi
1655 >           rij=xj*xj+yj*yj+zj*zj
1656 >           if (rij.lt.r0ijsq) then
1657 >             evdw1ij=0.25d0*(rij-r0ijsq)**2
1658 >             fac=rij-r0ijsq
1659 >           else
1660 >             evdw1ij=0.0d0
1661 >             fac=0.0d0
1662 1279,1288c1900,1906
1663 <       enddo
1664 <       if (calc_grad) then
1665 <       do i=1,nres-1
1666 <         vbld_inv_temp(1)=vbld_inv(i+1)
1667 <         if (i.lt.nres-1) then
1668 <           vbld_inv_temp(2)=vbld_inv(i+2)
1669 <         else
1670 <           vbld_inv_temp(2)=vbld_inv(i)
1671 <         endif
1672 <         do j=1,2
1673 ---
1674 >           evdw1=evdw1+evdw1ij
1675 > C
1676 > C Calculate contributions to the Cartesian gradient.
1677 > C
1678 >           ggg(1)=fac*xj
1679 >           ggg(2)=fac*yj
1680 >           ggg(3)=fac*zj
1681 1290,1293c1908,1909
1682 <             do l=1,3
1683 <               uygrad(l,k,j,i)=vbld_inv_temp(j)*uygrad(l,k,j,i)
1684 <               uzgrad(l,k,j,i)=vbld_inv_temp(j)*uzgrad(l,k,j,i)
1685 <             enddo
1686 ---
1687 >             gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1688 >             gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1689 1295,1297c1911,1930
1690 <         enddo
1691 <       enddo
1692 <       endif
1693 ---
1694 > *
1695 > * Loop over residues i+1 thru j-1.
1696 > *
1697 > cgrad          do k=i+1,j-1
1698 > cgrad            do l=1,3
1699 > cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
1700 > cgrad            enddo
1701 > cgrad          enddo
1702 >         enddo ! j
1703 >       enddo   ! i
1704 > cgrad      do i=nnt,nct-1
1705 > cgrad        do k=1,3
1706 > cgrad          gelc(k,i)=gelc(k,i)+0.5d0*gelc(k,i)
1707 > cgrad        enddo
1708 > cgrad        do j=i+1,nct-1
1709 > cgrad          do k=1,3
1710 > cgrad            gelc(k,i)=gelc(k,i)+gelc(k,j)
1711 > cgrad          enddo
1712 > cgrad        enddo
1713 > cgrad      enddo
1714 1300,1301c1933,1934
1715 < C-----------------------------------------------------------------------------
1716 <       subroutine vec_and_deriv_test
1717 ---
1718 > c------------------------------------------------------------------------------
1719 >       subroutine vec_and_deriv
1720 1304c1937,1939
1721 <       include 'DIMENSIONS.ZSCOPT'
1722 ---
1723 > #ifdef MPI
1724 >       include 'mpif.h'
1725 > #endif
1726 1311c1946,1948
1727 <       dimension uyder(3,3,2),uzder(3,3,2)
1728 ---
1729 >       include 'COMMON.SETUP'
1730 >       include 'COMMON.TIME1'
1731 >       dimension uyder(3,3,2),uzder(3,3,2),vbld_inv_temp(2)
1732 1314a1952,1954
1733 > #ifdef PARVEC
1734 >       do i=ivec_start,ivec_end
1735 > #else
1736 1315a1956
1737 > #endif
1738 1322,1324d1962
1739 < c            write (iout,*) 'fac',fac,
1740 < c     &        1.0d0/dsqrt(scalar(uz(1,i),uz(1,i)))
1741 <             fac=1.0d0/dsqrt(scalar(uz(1,i),uz(1,i)))
1742 1348,1350d1985
1743 <             do k=1,3
1744 <               uy(k,i)=fac*(dc_norm(k,i-1)-costh*dc_norm(k,i))
1745 <             enddo
1746 1352,1365d1986
1747 <             facy=1.0d0/dsqrt(scalar(dc_norm(1,i),dc_norm(1,i))*
1748 <      &       (scalar(dc_norm(1,i-1),dc_norm(1,i-1))**2-
1749 <      &        scalar(dc_norm(1,i),dc_norm(1,i-1))**2))
1750 <             do k=1,3
1751 < c              uy(k,i)=facy*(dc_norm(k,i+1)-costh*dc_norm(k,i))
1752 <               uy(k,i)=
1753 < c     &        facy*(
1754 <      &        dc_norm(k,i-1)*scalar(dc_norm(1,i),dc_norm(1,i))
1755 <      &        -scalar(dc_norm(1,i),dc_norm(1,i-1))*dc_norm(k,i)
1756 < c     &        )
1757 <             enddo
1758 < c            write (iout,*) 'facy',facy,
1759 < c     &       1.0d0/dsqrt(scalar(uy(1,i),uy(1,i)))
1760 <             facy=1.0d0/dsqrt(scalar(uy(1,i),uy(1,i)))
1761 1367c1988
1762 <               uy(k,i)=facy*uy(k,i)
1763 ---
1764 >               uy(k,i)=fac*(dc_norm(k,i-1)-costh*dc_norm(k,i))
1765 1376,1381c1997,1998
1766 < c              uyder(j,j,1)=uyder(j,j,1)-costh
1767 < c              uyder(j,j,2)=1.0d0+uyder(j,j,2)
1768 <               uyder(j,j,1)=uyder(j,j,1)
1769 <      &          -scalar(dc_norm(1,i),dc_norm(1,i-1))
1770 <               uyder(j,j,2)=scalar(dc_norm(1,i),dc_norm(1,i))
1771 <      &          +uyder(j,j,2)
1772 ---
1773 >               uyder(j,j,1)=uyder(j,j,1)-costh
1774 >               uyder(j,j,2)=1.0d0+uyder(j,j,2)
1775 1401d2017
1776 <             fac=1.0d0/dsqrt(scalar(uz(1,i),uz(1,i)))
1777 1426,1439d2041
1778 <             facy=1.0d0/dsqrt(scalar(dc_norm(1,i),dc_norm(1,i))*
1779 <      &       (scalar(dc_norm(1,i+1),dc_norm(1,i+1))**2-
1780 <      &        scalar(dc_norm(1,i),dc_norm(1,i+1))**2))
1781 <             do k=1,3
1782 < c              uy(k,i)=facy*(dc_norm(k,i+1)-costh*dc_norm(k,i))
1783 <               uy(k,i)=
1784 < c     &        facy*(
1785 <      &        dc_norm(k,i+1)*scalar(dc_norm(1,i),dc_norm(1,i))
1786 <      &        -scalar(dc_norm(1,i),dc_norm(1,i+1))*dc_norm(k,i)
1787 < c     &        )
1788 <             enddo
1789 < c            write (iout,*) 'facy',facy,
1790 < c     &       1.0d0/dsqrt(scalar(uy(1,i),uy(1,i)))
1791 <             facy=1.0d0/dsqrt(scalar(uy(1,i),uy(1,i)))
1792 1441c2043
1793 <               uy(k,i)=facy*uy(k,i)
1794 ---
1795 >               uy(k,i)=facy*(dc_norm(k,i+1)-costh*dc_norm(k,i))
1796 1450,1455c2052,2053
1797 < c              uyder(j,j,1)=uyder(j,j,1)-costh
1798 < c              uyder(j,j,2)=1.0d0+uyder(j,j,2)
1799 <               uyder(j,j,1)=uyder(j,j,1)
1800 <      &          -scalar(dc_norm(1,i),dc_norm(1,i+1))
1801 <               uyder(j,j,2)=scalar(dc_norm(1,i),dc_norm(1,i))
1802 <      &          +uyder(j,j,2)
1803 ---
1804 >               uyder(j,j,1)=uyder(j,j,1)-costh
1805 >               uyder(j,j,2)=1.0d0+uyder(j,j,2)
1806 1471a2070,2075
1807 >         vbld_inv_temp(1)=vbld_inv(i+1)
1808 >         if (i.lt.nres-1) then
1809 >           vbld_inv_temp(2)=vbld_inv(i+2)
1810 >           else
1811 >           vbld_inv_temp(2)=vbld_inv(i)
1812 >           endif
1813 1475,1476c2079,2080
1814 <               uygrad(l,k,j,i)=vblinv*uygrad(l,k,j,i)
1815 <               uzgrad(l,k,j,i)=vblinv*uzgrad(l,k,j,i)
1816 ---
1817 >               uygrad(l,k,j,i)=vbld_inv_temp(j)*uygrad(l,k,j,i)
1818 >               uzgrad(l,k,j,i)=vbld_inv_temp(j)*uzgrad(l,k,j,i)
1819 1480a2085,2112
1820 > #if defined(PARVEC) && defined(MPI)
1821 >       if (nfgtasks1.gt.1) then
1822 >         time00=MPI_Wtime()
1823 > c        print *,"Processor",fg_rank1,kolor1," ivec_start",ivec_start,
1824 > c     &   " ivec_displ",(ivec_displ(i),i=0,nfgtasks1-1),
1825 > c     &   " ivec_count",(ivec_count(i),i=0,nfgtasks1-1)
1826 >         call MPI_Allgatherv(uy(1,ivec_start),ivec_count(fg_rank1),
1827 >      &   MPI_UYZ,uy(1,1),ivec_count(0),ivec_displ(0),MPI_UYZ,
1828 >      &   FG_COMM1,IERR)
1829 >         call MPI_Allgatherv(uz(1,ivec_start),ivec_count(fg_rank1),
1830 >      &   MPI_UYZ,uz(1,1),ivec_count(0),ivec_displ(0),MPI_UYZ,
1831 >      &   FG_COMM1,IERR)
1832 >         call MPI_Allgatherv(uygrad(1,1,1,ivec_start),
1833 >      &   ivec_count(fg_rank1),MPI_UYZGRAD,uygrad(1,1,1,1),ivec_count(0),
1834 >      &   ivec_displ(0),MPI_UYZGRAD,FG_COMM1,IERR)
1835 >         call MPI_Allgatherv(uzgrad(1,1,1,ivec_start),
1836 >      &   ivec_count(fg_rank1),MPI_UYZGRAD,uzgrad(1,1,1,1),ivec_count(0),
1837 >      &   ivec_displ(0),MPI_UYZGRAD,FG_COMM1,IERR)
1838 >         time_gather=time_gather+MPI_Wtime()-time00
1839 >       endif
1840 > c      if (fg_rank.eq.0) then
1841 > c        write (iout,*) "Arrays UY and UZ"
1842 > c        do i=1,nres-1
1843 > c          write (iout,'(i5,3f10.5,5x,3f10.5)') i,(uy(k,i),k=1,3),
1844 > c     &     (uz(k,i),k=1,3)
1845 > c        enddo
1846 > c      endif
1847 > #endif
1848 1487d2118
1849 <       include 'DIMENSIONS.ZSCOPT'
1850 1572c2203,2208
1851 <       include 'DIMENSIONS.ZSCOPT'
1852 ---
1853 > #ifdef MPI
1854 >       include "mpif.h"
1855 >       include "COMMON.SETUP"
1856 >       integer IERR
1857 >       integer status(MPI_STATUS_SIZE)
1858 > #endif
1859 1588a2225,2227
1860 > #ifdef PARMAT
1861 >       do i=ivec_start+2,ivec_end+2
1862 > #else
1863 1589a2229
1864 > #endif
1865 1655a2296
1866 > c        if (i.gt. iatel_s+2 .and. i.lt.iatel_e+5) then
1867 1657,1661c2298
1868 <           if (itype(i-2).le.ntyp) then
1869 <             iti = itortyp(itype(i-2))
1870 <           else
1871 <             iti=ntortyp+1
1872 <           endif
1873 ---
1874 >           iti = itortyp(itype(i-2))
1875 1664a2302
1876 > c        if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
1877 1666,1670c2304
1878 <           if (itype(i-1).le.ntyp) then
1879 <             iti1 = itortyp(itype(i-1))
1880 <           else
1881 <             iti1=ntortyp+1
1882 <           endif
1883 ---
1884 >           iti1 = itortyp(itype(i-1))
1885 1678,1679c2312,2313
1886 < c        print *,"itilde1 i iti iti1",i,iti,iti1
1887 <         if (i .gt. iatel_s+2) then
1888 ---
1889 > c        if (i .gt. iatel_s+2) then
1890 >         if (i .gt. nnt+2) then
1891 1681a2316,2317
1892 >           if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0)
1893 >      &    then
1894 1686a2323
1895 >           endif
1896 1700d2336
1897 < c        print *,"itilde2 i iti iti1",i,iti,iti1
1898 1703,1708d2338
1899 <         call matmat2(CC(1,1,iti1),Ugder(1,1,i-2),CUgder(1,1,i-2))
1900 <         call matmat2(DD(1,1,iti),Ugder(1,1,i-2),DUgder(1,1,i-2))
1901 <         call matmat2(Dtilde(1,1,iti),Ug2der(1,1,i-2),DtUg2der(1,1,i-2))
1902 <         call matvec2(Ctilde(1,1,iti1),obrot_der(1,i-2),Ctobrder(1,i-2))
1903 <         call matvec2(Dtilde(1,1,iti),obrot2_der(1,i-2),Dtobr2der(1,i-2))
1904 < c        print *,"itilde3 i iti iti1",i,iti,iti1
1905 1711a2342
1906 > c        if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
1907 1713,1717c2344
1908 <           if (itype(i-1).le.ntyp) then
1909 <             iti1 = itortyp(itype(i-1))
1910 <           else
1911 <             iti1=ntortyp+1
1912 <           endif
1913 ---
1914 >           iti1 = itortyp(itype(i-1))
1915 1723a2351,2360
1916 > cd        write (iout,*) 'mu ',mu(:,i-2)
1917 > cd        write (iout,*) 'mu1',mu1(:,i-2)
1918 > cd        write (iout,*) 'mu2',mu2(:,i-2)
1919 >         if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0)
1920 >      &  then
1921 >         call matmat2(CC(1,1,iti1),Ugder(1,1,i-2),CUgder(1,1,i-2))
1922 >         call matmat2(DD(1,1,iti),Ugder(1,1,i-2),DUgder(1,1,i-2))
1923 >         call matmat2(Dtilde(1,1,iti),Ug2der(1,1,i-2),DtUg2der(1,1,i-2))
1924 >         call matvec2(Ctilde(1,1,iti1),obrot_der(1,i-2),Ctobrder(1,i-2))
1925 >         call matvec2(Dtilde(1,1,iti),obrot2_der(1,i-2),Dtobr2der(1,i-2))
1926 1734,1735c2371
1927 < cd        write (iout,*) 'i',i,' mu ',(mu(k,i-2),k=1,2),
1928 < cd     &  ' mu1',(b1(k,i-2),k=1,2),' mu2',(Ub2(k,i-2),k=1,2)
1929 ---
1930 >         endif
1931 1738a2375,2377
1932 >       if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0)
1933 >      &then
1934 > c      do i=max0(ivec_start,2),ivec_end
1935 1748a2388,2636
1936 >       endif
1937 > #if defined(MPI) && defined(PARMAT)
1938 > #ifdef DEBUG
1939 > c      if (fg_rank.eq.0) then
1940 >         write (iout,*) "Arrays UG and UGDER before GATHER"
1941 >         do i=1,nres-1
1942 >           write (iout,'(i5,4f10.5,5x,4f10.5)') i,
1943 >      &     ((ug(l,k,i),l=1,2),k=1,2),
1944 >      &     ((ugder(l,k,i),l=1,2),k=1,2)
1945 >         enddo
1946 >         write (iout,*) "Arrays UG2 and UG2DER"
1947 >         do i=1,nres-1
1948 >           write (iout,'(i5,4f10.5,5x,4f10.5)') i,
1949 >      &     ((ug2(l,k,i),l=1,2),k=1,2),
1950 >      &     ((ug2der(l,k,i),l=1,2),k=1,2)
1951 >         enddo
1952 >         write (iout,*) "Arrays OBROT OBROT2 OBROTDER and OBROT2DER"
1953 >         do i=1,nres-1
1954 >           write (iout,'(i5,4f10.5,5x,4f10.5)') i,
1955 >      &     (obrot(k,i),k=1,2),(obrot2(k,i),k=1,2),
1956 >      &     (obrot_der(k,i),k=1,2),(obrot2_der(k,i),k=1,2)
1957 >         enddo
1958 >         write (iout,*) "Arrays COSTAB SINTAB COSTAB2 and SINTAB2"
1959 >         do i=1,nres-1
1960 >           write (iout,'(i5,4f10.5,5x,4f10.5)') i,
1961 >      &     costab(i),sintab(i),costab2(i),sintab2(i)
1962 >         enddo
1963 >         write (iout,*) "Array MUDER"
1964 >         do i=1,nres-1
1965 >           write (iout,'(i5,2f10.5)') i,muder(1,i),muder(2,i)
1966 >         enddo
1967 > c      endif
1968 > #endif
1969 >       if (nfgtasks.gt.1) then
1970 >         time00=MPI_Wtime()
1971 > c        write(iout,*)"Processor",fg_rank,kolor," ivec_start",ivec_start,
1972 > c     &   " ivec_displ",(ivec_displ(i),i=0,nfgtasks-1),
1973 > c     &   " ivec_count",(ivec_count(i),i=0,nfgtasks-1)
1974 > #ifdef MATGATHER
1975 >         call MPI_Allgatherv(Ub2(1,ivec_start),ivec_count(fg_rank1),
1976 >      &   MPI_MU,Ub2(1,1),ivec_count(0),ivec_displ(0),MPI_MU,
1977 >      &   FG_COMM1,IERR)
1978 >         call MPI_Allgatherv(Ub2der(1,ivec_start),ivec_count(fg_rank1),
1979 >      &   MPI_MU,Ub2der(1,1),ivec_count(0),ivec_displ(0),MPI_MU,
1980 >      &   FG_COMM1,IERR)
1981 >         call MPI_Allgatherv(mu(1,ivec_start),ivec_count(fg_rank1),
1982 >      &   MPI_MU,mu(1,1),ivec_count(0),ivec_displ(0),MPI_MU,
1983 >      &   FG_COMM1,IERR)
1984 >         call MPI_Allgatherv(muder(1,ivec_start),ivec_count(fg_rank1),
1985 >      &   MPI_MU,muder(1,1),ivec_count(0),ivec_displ(0),MPI_MU,
1986 >      &   FG_COMM1,IERR)
1987 >         call MPI_Allgatherv(Eug(1,1,ivec_start),ivec_count(fg_rank1),
1988 >      &   MPI_MAT1,Eug(1,1,1),ivec_count(0),ivec_displ(0),MPI_MAT1,
1989 >      &   FG_COMM1,IERR)
1990 >         call MPI_Allgatherv(Eugder(1,1,ivec_start),ivec_count(fg_rank1),
1991 >      &   MPI_MAT1,Eugder(1,1,1),ivec_count(0),ivec_displ(0),MPI_MAT1,
1992 >      &   FG_COMM1,IERR)
1993 >         call MPI_Allgatherv(costab(ivec_start),ivec_count(fg_rank1),
1994 >      &   MPI_DOUBLE_PRECISION,costab(1),ivec_count(0),ivec_displ(0),
1995 >      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
1996 >         call MPI_Allgatherv(sintab(ivec_start),ivec_count(fg_rank1),
1997 >      &   MPI_DOUBLE_PRECISION,sintab(1),ivec_count(0),ivec_displ(0),
1998 >      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
1999 >         call MPI_Allgatherv(costab2(ivec_start),ivec_count(fg_rank1),
2000 >      &   MPI_DOUBLE_PRECISION,costab2(1),ivec_count(0),ivec_displ(0),
2001 >      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
2002 >         call MPI_Allgatherv(sintab2(ivec_start),ivec_count(fg_rank1),
2003 >      &   MPI_DOUBLE_PRECISION,sintab2(1),ivec_count(0),ivec_displ(0),
2004 >      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
2005 >         if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0)
2006 >      &  then
2007 >         call MPI_Allgatherv(Ctobr(1,ivec_start),ivec_count(fg_rank1),
2008 >      &   MPI_MU,Ctobr(1,1),ivec_count(0),ivec_displ(0),MPI_MU,
2009 >      &   FG_COMM1,IERR)
2010 >         call MPI_Allgatherv(Ctobrder(1,ivec_start),ivec_count(fg_rank1),
2011 >      &   MPI_MU,Ctobrder(1,1),ivec_count(0),ivec_displ(0),MPI_MU,
2012 >      &   FG_COMM1,IERR)
2013 >         call MPI_Allgatherv(Dtobr2(1,ivec_start),ivec_count(fg_rank1),
2014 >      &   MPI_MU,Dtobr2(1,1),ivec_count(0),ivec_displ(0),MPI_MU,
2015 >      &   FG_COMM1,IERR)
2016 >        call MPI_Allgatherv(Dtobr2der(1,ivec_start),ivec_count(fg_rank1),
2017 >      &   MPI_MU,Dtobr2der(1,1),ivec_count(0),ivec_displ(0),MPI_MU,
2018 >      &   FG_COMM1,IERR)
2019 >         call MPI_Allgatherv(Ug2Db1t(1,ivec_start),ivec_count(fg_rank1),
2020 >      &   MPI_MU,Ug2Db1t(1,1),ivec_count(0),ivec_displ(0),MPI_MU,
2021 >      &   FG_COMM1,IERR)
2022 >         call MPI_Allgatherv(Ug2Db1tder(1,ivec_start),
2023 >      &   ivec_count(fg_rank1),
2024 >      &   MPI_MU,Ug2Db1tder(1,1),ivec_count(0),ivec_displ(0),MPI_MU,
2025 >      &   FG_COMM1,IERR)
2026 >         call MPI_Allgatherv(CUgb2(1,ivec_start),ivec_count(fg_rank1),
2027 >      &   MPI_MU,CUgb2(1,1),ivec_count(0),ivec_displ(0),MPI_MU,
2028 >      &   FG_COMM1,IERR)
2029 >         call MPI_Allgatherv(CUgb2der(1,ivec_start),ivec_count(fg_rank1),
2030 >      &   MPI_MU,CUgb2der(1,1),ivec_count(0),ivec_displ(0),MPI_MU,
2031 >      &   FG_COMM1,IERR)
2032 >         call MPI_Allgatherv(Cug(1,1,ivec_start),ivec_count(fg_rank1),
2033 >      &   MPI_MAT1,Cug(1,1,1),ivec_count(0),ivec_displ(0),MPI_MAT1,
2034 >      &   FG_COMM1,IERR)
2035 >         call MPI_Allgatherv(Cugder(1,1,ivec_start),ivec_count(fg_rank1),
2036 >      &   MPI_MAT1,Cugder(1,1,1),ivec_count(0),ivec_displ(0),MPI_MAT1,
2037 >      &   FG_COMM1,IERR)
2038 >         call MPI_Allgatherv(Dug(1,1,ivec_start),ivec_count(fg_rank1),
2039 >      &   MPI_MAT1,Dug(1,1,1),ivec_count(0),ivec_displ(0),MPI_MAT1,
2040 >      &   FG_COMM1,IERR)
2041 >         call MPI_Allgatherv(Dugder(1,1,ivec_start),ivec_count(fg_rank1),
2042 >      &   MPI_MAT1,Dugder(1,1,1),ivec_count(0),ivec_displ(0),MPI_MAT1,
2043 >      &   FG_COMM1,IERR)
2044 >         call MPI_Allgatherv(Dtug2(1,1,ivec_start),ivec_count(fg_rank1),
2045 >      &   MPI_MAT1,Dtug2(1,1,1),ivec_count(0),ivec_displ(0),MPI_MAT1,
2046 >      &   FG_COMM1,IERR)
2047 >         call MPI_Allgatherv(Dtug2der(1,1,ivec_start),
2048 >      &   ivec_count(fg_rank1),
2049 >      &   MPI_MAT1,Dtug2der(1,1,1),ivec_count(0),ivec_displ(0),MPI_MAT1,
2050 >      &   FG_COMM1,IERR)
2051 >         call MPI_Allgatherv(EugC(1,1,ivec_start),ivec_count(fg_rank1),
2052 >      &   MPI_MAT1,EugC(1,1,1),ivec_count(0),ivec_displ(0),MPI_MAT1,
2053 >      &   FG_COMM1,IERR)
2054 >        call MPI_Allgatherv(EugCder(1,1,ivec_start),ivec_count(fg_rank1),
2055 >      &   MPI_MAT1,EugCder(1,1,1),ivec_count(0),ivec_displ(0),MPI_MAT1,
2056 >      &   FG_COMM1,IERR)
2057 >         call MPI_Allgatherv(EugD(1,1,ivec_start),ivec_count(fg_rank1),
2058 >      &   MPI_MAT1,EugD(1,1,1),ivec_count(0),ivec_displ(0),MPI_MAT1,
2059 >      &   FG_COMM1,IERR)
2060 >        call MPI_Allgatherv(EugDder(1,1,ivec_start),ivec_count(fg_rank1),
2061 >      &   MPI_MAT1,EugDder(1,1,1),ivec_count(0),ivec_displ(0),MPI_MAT1,
2062 >      &   FG_COMM1,IERR)
2063 >         call MPI_Allgatherv(DtUg2EUg(1,1,ivec_start),
2064 >      &   ivec_count(fg_rank1),
2065 >      &   MPI_MAT1,DtUg2EUg(1,1,1),ivec_count(0),ivec_displ(0),MPI_MAT1,
2066 >      &   FG_COMM1,IERR)
2067 >         call MPI_Allgatherv(Ug2DtEUg(1,1,ivec_start),
2068 >      &   ivec_count(fg_rank1),
2069 >      &   MPI_MAT1,Ug2DtEUg(1,1,1),ivec_count(0),ivec_displ(0),MPI_MAT1,
2070 >      &   FG_COMM1,IERR)
2071 >         call MPI_Allgatherv(DtUg2EUgder(1,1,1,ivec_start),
2072 >      &   ivec_count(fg_rank1),
2073 >      &   MPI_MAT2,DtUg2EUgder(1,1,1,1),ivec_count(0),ivec_displ(0),
2074 >      &   MPI_MAT2,FG_COMM1,IERR)
2075 >         call MPI_Allgatherv(Ug2DtEUgder(1,1,1,ivec_start),
2076 >      &   ivec_count(fg_rank1),
2077 >      &   MPI_MAT2,Ug2DtEUgder(1,1,1,1),ivec_count(0),ivec_displ(0),
2078 >      &   MPI_MAT2,FG_COMM1,IERR)
2079 >         endif
2080 > #else
2081 > c Passes matrix info through the ring
2082 >       isend=fg_rank1
2083 >       irecv=fg_rank1-1
2084 >       if (irecv.lt.0) irecv=nfgtasks1-1
2085 >       iprev=irecv
2086 >       inext=fg_rank1+1
2087 >       if (inext.ge.nfgtasks1) inext=0
2088 >       do i=1,nfgtasks1-1
2089 > c        write (iout,*) "isend",isend," irecv",irecv
2090 > c        call flush(iout)
2091 >         lensend=lentyp(isend)
2092 >         lenrecv=lentyp(irecv)
2093 > c        write (iout,*) "lensend",lensend," lenrecv",lenrecv
2094 > c        call MPI_SENDRECV(ug(1,1,ivec_displ(isend)+1),1,
2095 > c     &   MPI_ROTAT1(lensend),inext,2200+isend,
2096 > c     &   ug(1,1,ivec_displ(irecv)+1),1,MPI_ROTAT1(lenrecv),
2097 > c     &   iprev,2200+irecv,FG_COMM,status,IERR)
2098 > c        write (iout,*) "Gather ROTAT1"
2099 > c        call flush(iout)
2100 > c        call MPI_SENDRECV(obrot(1,ivec_displ(isend)+1),1,
2101 > c     &   MPI_ROTAT2(lensend),inext,3300+isend,
2102 > c     &   obrot(1,ivec_displ(irecv)+1),1,MPI_ROTAT2(lenrecv),
2103 > c     &   iprev,3300+irecv,FG_COMM,status,IERR)
2104 > c        write (iout,*) "Gather ROTAT2"
2105 > c        call flush(iout)
2106 >         call MPI_SENDRECV(costab(ivec_displ(isend)+1),1,
2107 >      &   MPI_ROTAT_OLD(lensend),inext,4400+isend,
2108 >      &   costab(ivec_displ(irecv)+1),1,MPI_ROTAT_OLD(lenrecv),
2109 >      &   iprev,4400+irecv,FG_COMM,status,IERR)
2110 > c        write (iout,*) "Gather ROTAT_OLD"
2111 > c        call flush(iout)
2112 >         call MPI_SENDRECV(mu(1,ivec_displ(isend)+1),1,
2113 >      &   MPI_PRECOMP11(lensend),inext,5500+isend,
2114 >      &   mu(1,ivec_displ(irecv)+1),1,MPI_PRECOMP11(lenrecv),
2115 >      &   iprev,5500+irecv,FG_COMM,status,IERR)
2116 > c        write (iout,*) "Gather PRECOMP11"
2117 > c        call flush(iout)
2118 >         call MPI_SENDRECV(Eug(1,1,ivec_displ(isend)+1),1,
2119 >      &   MPI_PRECOMP12(lensend),inext,6600+isend,
2120 >      &   Eug(1,1,ivec_displ(irecv)+1),1,MPI_PRECOMP12(lenrecv),
2121 >      &   iprev,6600+irecv,FG_COMM,status,IERR)
2122 > c        write (iout,*) "Gather PRECOMP12"
2123 > c        call flush(iout)
2124 >         if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0)
2125 >      &  then
2126 >         call MPI_SENDRECV(ug2db1t(1,ivec_displ(isend)+1),1,
2127 >      &   MPI_ROTAT2(lensend),inext,7700+isend,
2128 >      &   ug2db1t(1,ivec_displ(irecv)+1),1,MPI_ROTAT2(lenrecv),
2129 >      &   iprev,7700+irecv,FG_COMM,status,IERR)
2130 > c        write (iout,*) "Gather PRECOMP21"
2131 > c        call flush(iout)
2132 >         call MPI_SENDRECV(EUgC(1,1,ivec_displ(isend)+1),1,
2133 >      &   MPI_PRECOMP22(lensend),inext,8800+isend,
2134 >      &   EUgC(1,1,ivec_displ(irecv)+1),1,MPI_PRECOMP22(lenrecv),
2135 >      &   iprev,8800+irecv,FG_COMM,status,IERR)
2136 > c        write (iout,*) "Gather PRECOMP22"
2137 > c        call flush(iout)
2138 >         call MPI_SENDRECV(Ug2DtEUgder(1,1,1,ivec_displ(isend)+1),1,
2139 >      &   MPI_PRECOMP23(lensend),inext,9900+isend,
2140 >      &   Ug2DtEUgder(1,1,1,ivec_displ(irecv)+1),1,
2141 >      &   MPI_PRECOMP23(lenrecv),
2142 >      &   iprev,9900+irecv,FG_COMM,status,IERR)
2143 > c        write (iout,*) "Gather PRECOMP23"
2144 > c        call flush(iout)
2145 >         endif
2146 >         isend=irecv
2147 >         irecv=irecv-1
2148 >         if (irecv.lt.0) irecv=nfgtasks1-1
2149 >       enddo
2150 > #endif
2151 >         time_gather=time_gather+MPI_Wtime()-time00
2152 >       endif
2153 > #ifdef DEBUG
2154 > c      if (fg_rank.eq.0) then
2155 >         write (iout,*) "Arrays UG and UGDER"
2156 >         do i=1,nres-1
2157 >           write (iout,'(i5,4f10.5,5x,4f10.5)') i,
2158 >      &     ((ug(l,k,i),l=1,2),k=1,2),
2159 >      &     ((ugder(l,k,i),l=1,2),k=1,2)
2160 >         enddo
2161 >         write (iout,*) "Arrays UG2 and UG2DER"
2162 >         do i=1,nres-1
2163 >           write (iout,'(i5,4f10.5,5x,4f10.5)') i,
2164 >      &     ((ug2(l,k,i),l=1,2),k=1,2),
2165 >      &     ((ug2der(l,k,i),l=1,2),k=1,2)
2166 >         enddo
2167 >         write (iout,*) "Arrays OBROT OBROT2 OBROTDER and OBROT2DER"
2168 >         do i=1,nres-1
2169 >           write (iout,'(i5,4f10.5,5x,4f10.5)') i,
2170 >      &     (obrot(k,i),k=1,2),(obrot2(k,i),k=1,2),
2171 >      &     (obrot_der(k,i),k=1,2),(obrot2_der(k,i),k=1,2)
2172 >         enddo
2173 >         write (iout,*) "Arrays COSTAB SINTAB COSTAB2 and SINTAB2"
2174 >         do i=1,nres-1
2175 >           write (iout,'(i5,4f10.5,5x,4f10.5)') i,
2176 >      &     costab(i),sintab(i),costab2(i),sintab2(i)
2177 >         enddo
2178 >         write (iout,*) "Array MUDER"
2179 >         do i=1,nres-1
2180 >           write (iout,'(i5,2f10.5)') i,muder(1,i),muder(2,i)
2181 >         enddo
2182 > c      endif
2183 > #endif
2184 > #endif
2185 1768a2657,2659
2186 > #ifdef MPI
2187 >       include 'mpif.h'
2188 > #endif
2189 1770d2660
2190 <       include 'DIMENSIONS.ZSCOPT'
2191 1771a2662
2192 >       include 'COMMON.SETUP'
2193 1782a2674
2194 >       include 'COMMON.TIME1'
2195 1787c2679,2681
2196 <       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,j1
2197 ---
2198 >       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
2199 >      &    dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
2200 >      &    num_conti,j1,j2
2201 1788a2683,2685
2202 > #ifdef MOMENT
2203 >       double precision scal_el /1.0d0/
2204 > #else
2205 1789a2687
2206 > #endif
2207 1818,1823c2716,2719
2208 < cd      if (wel_loc.gt.0.0d0) then
2209 <         if (icheckgrad.eq.1) then
2210 <         call vec_and_deriv_test
2211 <         else
2212 <         call vec_and_deriv
2213 <         endif
2214 ---
2215 > c        call vec_and_deriv
2216 > #ifdef TIMING
2217 >         time01=MPI_Wtime()
2218 > #endif
2219 1824a2721,2723
2220 > #ifdef TIMING
2221 >         time_mat=time_mat+MPI_Wtime()-time01
2222 > #endif
2223 1829c2728
2224 < cd          write (iout,'(i5,2f10.5)') k,uy(k,i),uz(k,i)
2225 ---
2226 > cd        write (iout,'(i5,2f10.5)') k,uy(k,i),uz(k,i)
2227 1836c2735
2228 <       num_conti_hb=0
2229 ---
2230 >       t_eelecij=0.0d0
2231 1851a2751,2795
2232 > c
2233 > c
2234 > c 9/27/08 AL Split the interaction loop to ensure load balancing of turn terms
2235 > C
2236 > C Loop over i,i+2 and i,i+3 pairs of the peptide groups
2237 > C
2238 >       do i=iturn3_start,iturn3_end
2239 >         if (itype(i).eq.21 .or. itype(i+1).eq.21
2240 >      &  .or. itype(i+2).eq.21 .or. itype(i+3).eq.21) cycle
2241 >         dxi=dc(1,i)
2242 >         dyi=dc(2,i)
2243 >         dzi=dc(3,i)
2244 >         dx_normi=dc_norm(1,i)
2245 >         dy_normi=dc_norm(2,i)
2246 >         dz_normi=dc_norm(3,i)
2247 >         xmedi=c(1,i)+0.5d0*dxi
2248 >         ymedi=c(2,i)+0.5d0*dyi
2249 >         zmedi=c(3,i)+0.5d0*dzi
2250 >         num_conti=0
2251 >         call eelecij(i,i+2,ees,evdw1,eel_loc)
2252 >         if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
2253 >         num_cont_hb(i)=num_conti
2254 >       enddo
2255 >       do i=iturn4_start,iturn4_end
2256 >         if (itype(i).eq.21 .or. itype(i+1).eq.21
2257 >      &    .or. itype(i+3).eq.21
2258 >      &    .or. itype(i+4).eq.21) cycle
2259 >         dxi=dc(1,i)
2260 >         dyi=dc(2,i)
2261 >         dzi=dc(3,i)
2262 >         dx_normi=dc_norm(1,i)
2263 >         dy_normi=dc_norm(2,i)
2264 >         dz_normi=dc_norm(3,i)
2265 >         xmedi=c(1,i)+0.5d0*dxi
2266 >         ymedi=c(2,i)+0.5d0*dyi
2267 >         zmedi=c(3,i)+0.5d0*dzi
2268 >         num_conti=num_cont_hb(i)
2269 >         call eelecij(i,i+3,ees,evdw1,eel_loc)
2270 >         if (wturn4.gt.0.0d0 .and. itype(i+2).ne.21)
2271 >      &   call eturn4(i,eello_turn4)
2272 >         num_cont_hb(i)=num_conti
2273 >       enddo   ! i
2274 > c
2275 > c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
2276 > c
2277 1854d2797
2278 <         if (itel(i).eq.0) goto 1215
2279 1864d2806
2280 <         num_conti=0
2281 1865a2808
2282 >         num_conti=num_cont_hb(i)
2283 1866a2810
2284 > c          write (iout,*) i,j,itype(i),itype(j)
2285 1868,1869c2812,2866
2286 <           if (itel(j).eq.0) goto 1216
2287 <           ind=ind+1
2288 ---
2289 >           call eelecij(i,j,ees,evdw1,eel_loc)
2290 >         enddo ! j
2291 >         num_cont_hb(i)=num_conti
2292 >       enddo   ! i
2293 > c      write (iout,*) "Number of loop steps in EELEC:",ind
2294 > cd      do i=1,nres
2295 > cd        write (iout,'(i3,3f10.5,5x,3f10.5)')
2296 > cd     &     i,(gel_loc(k,i),k=1,3),gel_loc_loc(i)
2297 > cd      enddo
2298 > c 12/7/99 Adam eello_turn3 will be considered as a separate energy term
2299 > ccc      eel_loc=eel_loc+eello_turn3
2300 > cd      print *,"Processor",fg_rank," t_eelecij",t_eelecij
2301 >       return
2302 >       end
2303 > C-------------------------------------------------------------------------------
2304 >       subroutine eelecij(i,j,ees,evdw1,eel_loc)
2305 >       implicit real*8 (a-h,o-z)
2306 >       include 'DIMENSIONS'
2307 > #ifdef MPI
2308 >       include "mpif.h"
2309 > #endif
2310 >       include 'COMMON.CONTROL'
2311 >       include 'COMMON.IOUNITS'
2312 >       include 'COMMON.GEO'
2313 >       include 'COMMON.VAR'
2314 >       include 'COMMON.LOCAL'
2315 >       include 'COMMON.CHAIN'
2316 >       include 'COMMON.DERIV'
2317 >       include 'COMMON.INTERACT'
2318 >       include 'COMMON.CONTACTS'
2319 >       include 'COMMON.TORSION'
2320 >       include 'COMMON.VECTORS'
2321 >       include 'COMMON.FFIELD'
2322 >       include 'COMMON.TIME1'
2323 >       dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
2324 >      &          erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
2325 >       double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
2326 >      &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
2327 >       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
2328 >      &    dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
2329 >      &    num_conti,j1,j2
2330 > c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
2331 > #ifdef MOMENT
2332 >       double precision scal_el /1.0d0/
2333 > #else
2334 >       double precision scal_el /0.5d0/
2335 > #endif
2336 > C 12/13/98
2337 > C 13-go grudnia roku pamietnego...
2338 >       double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
2339 >      &                   0.0d0,1.0d0,0.0d0,
2340 >      &                   0.0d0,0.0d0,1.0d0/
2341 > c          time00=MPI_Wtime()
2342 > cd      write (iout,*) "eelecij",i,j
2343 > c          ind=ind+1
2344 1875,1880d2871
2345 < C Diagnostics only!!!
2346 < c         aaa=0.0D0
2347 < c         bbb=0.0D0
2348 < c         ael6i=0.0D0
2349 < c         ael3i=0.0D0
2350 < C End diagnostics
2351 1912d2902
2352 < c          write (iout,*) "i",i,iteli," j",j,itelj," eesij",eesij
2353 1920a2911,2916
2354 >
2355 >           if (energy_dec) then
2356 >               write (iout,'(a6,2i5,0pf7.3)') 'evdw1',i,j,evdwij
2357 >               write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
2358 >           endif
2359 >
2360 1925c2921
2361 <           facvdw=-6*rrmij*(ev1+evdwij)
2362 ---
2363 >           facvdw=-6*rrmij*(ev1+evdwij)
2364 1931d2926
2365 <           if (calc_grad) then
2366 1934c2929
2367 < *
2368 ---
2369 > *
2370 1937a2933,2938
2371 > c          do k=1,3
2372 > c            ghalf=0.5D0*ggg(k)
2373 > c            gelc(k,i)=gelc(k,i)+ghalf
2374 > c            gelc(k,j)=gelc(k,j)+ghalf
2375 > c          enddo
2376 > c 9/28/08 AL Gradient compotents will be summed only at the end
2377 1939,1941c2940,2941
2378 <             ghalf=0.5D0*ggg(k)
2379 <             gelc(k,i)=gelc(k,i)+ghalf
2380 <             gelc(k,j)=gelc(k,j)+ghalf
2381 ---
2382 >             gelc_long(k,j)=gelc_long(k,j)+ggg(k)
2383 >             gelc_long(k,i)=gelc_long(k,i)-ggg(k)
2384 1946,1950c2946,2950
2385 <           do k=i+1,j-1
2386 <             do l=1,3
2387 <               gelc(l,k)=gelc(l,k)+ggg(l)
2388 <             enddo
2389 <           enddo
2390 ---
2391 > cgrad          do k=i+1,j-1
2392 > cgrad            do l=1,3
2393 > cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
2394 > cgrad            enddo
2395 > cgrad          enddo
2396 1953a2954,2959
2397 > c          do k=1,3
2398 > c            ghalf=0.5D0*ggg(k)
2399 > c            gvdwpp(k,i)=gvdwpp(k,i)+ghalf
2400 > c            gvdwpp(k,j)=gvdwpp(k,j)+ghalf
2401 > c          enddo
2402 > c 9/28/08 AL Gradient compotents will be summed only at the end
2403 1955,1957c2961,2962
2404 <             ghalf=0.5D0*ggg(k)
2405 <             gvdwpp(k,i)=gvdwpp(k,i)+ghalf
2406 <             gvdwpp(k,j)=gvdwpp(k,j)+ghalf
2407 ---
2408 >             gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
2409 >             gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
2410 1962,1966c2967,2971
2411 <           do k=i+1,j-1
2412 <             do l=1,3
2413 <               gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
2414 <             enddo
2415 <           enddo
2416 ---
2417 > cgrad          do k=i+1,j-1
2418 > cgrad            do l=1,3
2419 > cgrad              gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
2420 > cgrad            enddo
2421 > cgrad          enddo
2422 1975d2979
2423 <           if (calc_grad) then
2424 1981a2986,2991
2425 > c          do k=1,3
2426 > c            ghalf=0.5D0*ggg(k)
2427 > c            gelc(k,i)=gelc(k,i)+ghalf
2428 > c            gelc(k,j)=gelc(k,j)+ghalf
2429 > c          enddo
2430 > c 9/28/08 AL Gradient compotents will be summed only at the end
2431 1983,1985c2993,2994
2432 <             ghalf=0.5D0*ggg(k)
2433 <             gelc(k,i)=gelc(k,i)+ghalf
2434 <             gelc(k,j)=gelc(k,j)+ghalf
2435 ---
2436 >             gelc_long(k,j)=gelc(k,j)+ggg(k)
2437 >             gelc_long(k,i)=gelc(k,i)-ggg(k)
2438 1990,1993c2999,3010
2439 <           do k=i+1,j-1
2440 <             do l=1,3
2441 <               gelc(l,k)=gelc(l,k)+ggg(l)
2442 <             enddo
2443 ---
2444 > cgrad          do k=i+1,j-1
2445 > cgrad            do l=1,3
2446 > cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
2447 > cgrad            enddo
2448 > cgrad          enddo
2449 > c 9/28/08 AL Gradient compotents will be summed only at the end
2450 >           ggg(1)=facvdw*xj
2451 >           ggg(2)=facvdw*yj
2452 >           ggg(3)=facvdw*zj
2453 >           do k=1,3
2454 >             gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
2455 >             gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
2456 2012a3030,3043
2457 > c          do k=1,3
2458 > c            ghalf=0.5D0*ggg(k)
2459 > c            gelc(k,i)=gelc(k,i)+ghalf
2460 > c     &               +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
2461 > c     &               + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2462 > c            gelc(k,j)=gelc(k,j)+ghalf
2463 > c     &               +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
2464 > c     &               + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2465 > c          enddo
2466 > cgrad          do k=i+1,j-1
2467 > cgrad            do l=1,3
2468 > cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
2469 > cgrad            enddo
2470 > cgrad          enddo
2471 2014,2015c3045
2472 <             ghalf=0.5D0*ggg(k)
2473 <             gelc(k,i)=gelc(k,i)+ghalf
2474 ---
2475 >             gelc(k,i)=gelc(k,i)
2476 2018c3048
2477 <             gelc(k,j)=gelc(k,j)+ghalf
2478 ---
2479 >             gelc(k,j)=gelc(k,j)
2480 2020a3051,3052
2481 >             gelc_long(k,j)=gelc_long(k,j)+ggg(k)
2482 >             gelc_long(k,i)=gelc_long(k,i)-ggg(k)
2483 2022,2028d3053
2484 <           do k=i+1,j-1
2485 <             do l=1,3
2486 <               gelc(l,k)=gelc(l,k)+ggg(l)
2487 <             enddo
2488 <           enddo
2489 <           endif
2490 <
2491 2064,2068d3088
2492 < C For diagnostics only
2493 < cd          a22=1.0d0
2494 < cd          a23=1.0d0
2495 < cd          a32=1.0d0
2496 < cd          a33=1.0d0
2497 2070,2072d3089
2498 < cd          write (2,*) 'fac=',fac
2499 < C For diagnostics only
2500 < cd          fac=1.0d0
2501 2080,2081c3097,3098
2502 < cd          write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') (uy(k,i),k=1,3),
2503 < cd     &      (uz(k,i),k=1,3),(uy(k,j),k=1,3),(uz(k,j),k=1,3)
2504 ---
2505 > cd          write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
2506 > cd     &      uy(:,j),uz(:,j)
2507 2086c3103
2508 < cd           write (iout,'(2i3,9f10.5/)') i,j,
2509 ---
2510 > cd           write (iout,'(9f10.5/)')
2511 2088d3104
2512 <           if (calc_grad) then
2513 2090,2095c3106
2514 <           call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
2515 < cd          do k=1,3
2516 < cd            do l=1,3
2517 < cd              erder(k,l)=0.0d0
2518 < cd            enddo
2519 < cd          enddo
2520 ---
2521 >           call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
2522 2110,2117d3120
2523 < cd          do k=1,3
2524 < cd            do l=1,3
2525 < cd              uryg(k,l)=0.0d0
2526 < cd              urzg(k,l)=0.0d0
2527 < cd              vryg(k,l)=0.0d0
2528 < cd              vrzg(k,l)=0.0d0
2529 < cd            enddo
2530 < cd          enddo
2531 2124,2127d3126
2532 < cd          a22der=0.0d0
2533 < cd          a23der=0.0d0
2534 < cd          a32der=0.0d0
2535 < cd          a33der=0.0d0
2536 2150,2153c3149,3152
2537 <             ghalf1=0.5d0*agg(k,1)
2538 <             ghalf2=0.5d0*agg(k,2)
2539 <             ghalf3=0.5d0*agg(k,3)
2540 <             ghalf4=0.5d0*agg(k,4)
2541 ---
2542 > cgrad            ghalf1=0.5d0*agg(k,1)
2543 > cgrad            ghalf2=0.5d0*agg(k,2)
2544 > cgrad            ghalf3=0.5d0*agg(k,3)
2545 > cgrad            ghalf4=0.5d0*agg(k,4)
2546 2155c3154
2547 <      &      -3.0d0*uryg(k,2)*vry)+ghalf1
2548 ---
2549 >      &      -3.0d0*uryg(k,2)*vry)!+ghalf1
2550 2157c3156
2551 <      &      -3.0d0*uryg(k,2)*vrz)+ghalf2
2552 ---
2553 >      &      -3.0d0*uryg(k,2)*vrz)!+ghalf2
2554 2159c3158
2555 <      &      -3.0d0*urzg(k,2)*vry)+ghalf3
2556 ---
2557 >      &      -3.0d0*urzg(k,2)*vry)!+ghalf3
2558 2161c3160
2559 <      &      -3.0d0*urzg(k,2)*vrz)+ghalf4
2560 ---
2561 >      &      -3.0d0*urzg(k,2)*vrz)!+ghalf4
2562 2164c3163
2563 <      &      -3.0d0*uryg(k,3)*vry)+agg(k,1)
2564 ---
2565 >      &      -3.0d0*uryg(k,3)*vry)!+agg(k,1)
2566 2166c3165
2567 <      &      -3.0d0*uryg(k,3)*vrz)+agg(k,2)
2568 ---
2569 >      &      -3.0d0*uryg(k,3)*vrz)!+agg(k,2)
2570 2168c3167
2571 <      &      -3.0d0*urzg(k,3)*vry)+agg(k,3)
2572 ---
2573 >      &      -3.0d0*urzg(k,3)*vry)!+agg(k,3)
2574 2170c3169
2575 <      &      -3.0d0*urzg(k,3)*vrz)+agg(k,4)
2576 ---
2577 >      &      -3.0d0*urzg(k,3)*vrz)!+agg(k,4)
2578 2173c3172
2579 <      &      -3.0d0*vryg(k,2)*ury)+ghalf1
2580 ---
2581 >      &      -3.0d0*vryg(k,2)*ury)!+ghalf1
2582 2175c3174
2583 <      &      -3.0d0*vrzg(k,2)*ury)+ghalf2
2584 ---
2585 >      &      -3.0d0*vrzg(k,2)*ury)!+ghalf2
2586 2177c3176
2587 <      &      -3.0d0*vryg(k,2)*urz)+ghalf3
2588 ---
2589 >      &      -3.0d0*vryg(k,2)*urz)!+ghalf3
2590 2179c3178
2591 <      &      -3.0d0*vrzg(k,2)*urz)+ghalf4
2592 ---
2593 >      &      -3.0d0*vrzg(k,2)*urz)!+ghalf4
2594 2189,2213c3188,3192
2595 < cd            aggi(k,1)=ghalf1
2596 < cd            aggi(k,2)=ghalf2
2597 < cd            aggi(k,3)=ghalf3
2598 < cd            aggi(k,4)=ghalf4
2599 < C Derivatives in DC(i+1)
2600 < cd            aggi1(k,1)=agg(k,1)
2601 < cd            aggi1(k,2)=agg(k,2)
2602 < cd            aggi1(k,3)=agg(k,3)
2603 < cd            aggi1(k,4)=agg(k,4)
2604 < C Derivatives in DC(j)
2605 < cd            aggj(k,1)=ghalf1
2606 < cd            aggj(k,2)=ghalf2
2607 < cd            aggj(k,3)=ghalf3
2608 < cd            aggj(k,4)=ghalf4
2609 < C Derivatives in DC(j+1)
2610 < cd            aggj1(k,1)=0.0d0
2611 < cd            aggj1(k,2)=0.0d0
2612 < cd            aggj1(k,3)=0.0d0
2613 < cd            aggj1(k,4)=0.0d0
2614 <             if (j.eq.nres-1 .and. i.lt.j-2) then
2615 <               do l=1,4
2616 <                 aggj1(k,l)=aggj1(k,l)+agg(k,l)
2617 < cd                aggj1(k,l)=agg(k,l)
2618 <               enddo
2619 <             endif
2620 ---
2621 > cgrad            if (j.eq.nres-1 .and. i.lt.j-2) then
2622 > cgrad              do l=1,4
2623 > cgrad                aggj1(k,l)=aggj1(k,l)+agg(k,l)
2624 > cgrad              enddo
2625 > cgrad            endif
2626 2215,2217d3193
2627 <           endif
2628 < c          goto 11111
2629 < C Check the loc-el terms by numerical integration
2630 2261d3236
2631 < 11111     continue
2632 2267c3242,3245
2633 < cd          write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3)
2634 ---
2635 >
2636 >           if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
2637 >      &            'eelloc',i,j,eel_loc_ij
2638 >
2639 2270d3247
2640 <           if (calc_grad) then
2641 2278,2284d3254
2642 < cd          call checkint3(i,j,mu1,mu2,a22,a23,a32,a33,acipa,eel_loc_ij)
2643 < cd          write(iout,*) 'agg  ',agg
2644 < cd          write(iout,*) 'aggi ',aggi
2645 < cd          write(iout,*) 'aggi1',aggi1
2646 < cd          write(iout,*) 'aggj ',aggj
2647 < cd          write(iout,*) 'aggj1',aggj1
2648 <
2649 2289,2294c3259,3269
2650 <           enddo
2651 <           do k=i+2,j2
2652 <             do l=1,3
2653 <               gel_loc(l,k)=gel_loc(l,k)+ggg(l)
2654 <             enddo
2655 <           enddo
2656 ---
2657 >             gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
2658 >             gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
2659 > cgrad            ghalf=0.5d0*ggg(l)
2660 > cgrad            gel_loc(l,i)=gel_loc(l,i)+ghalf
2661 > cgrad            gel_loc(l,j)=gel_loc(l,j)+ghalf
2662 >           enddo
2663 > cgrad          do k=i+1,j2
2664 > cgrad            do l=1,3
2665 > cgrad              gel_loc(l,k)=gel_loc(l,k)+ggg(l)
2666 > cgrad            enddo
2667 > cgrad          enddo
2668 2306d3280
2669 <           endif
2670 2308,2315d3281
2671 <           if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
2672 < C Contributions from turns
2673 <             a_temp(1,1)=a22
2674 <             a_temp(1,2)=a23
2675 <             a_temp(2,1)=a32
2676 <             a_temp(2,2)=a33
2677 <             call eturn34(i,j,eello_turn3,eello_turn4)
2678 <           endif
2679 2317c3283,3286
2680 <           if (j.gt.i+1 .and. num_conti.le.maxconts) then
2681 ---
2682 > c          if (j.gt.i+1 .and. num_conti.le.maxconts) then
2683 >           if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0
2684 >      &       .and. num_conti.le.maxconts) then
2685 > c            write (iout,*) i,j," entered corr"
2686 2334a3304,3305
2687 > cd                write (iout,*) "i",i," j",j," num_conti",num_conti,
2688 > cd     &           " jcont_hb",jcont_hb(num_conti,i)
2689 2350,2382d3320
2690 < c             if (i.eq.1) then
2691 < c                a_chuj(1,1,num_conti,i)=-0.61d0
2692 < c                a_chuj(1,2,num_conti,i)= 0.4d0
2693 < c                a_chuj(2,1,num_conti,i)= 0.65d0
2694 < c                a_chuj(2,2,num_conti,i)= 0.50d0
2695 < c             else if (i.eq.2) then
2696 < c                a_chuj(1,1,num_conti,i)= 0.0d0
2697 < c                a_chuj(1,2,num_conti,i)= 0.0d0
2698 < c                a_chuj(2,1,num_conti,i)= 0.0d0
2699 < c                a_chuj(2,2,num_conti,i)= 0.0d0
2700 < c             endif
2701 < C     --- and its gradients
2702 < cd                write (iout,*) 'i',i,' j',j
2703 < cd                do kkk=1,3
2704 < cd                write (iout,*) 'iii 1 kkk',kkk
2705 < cd                write (iout,*) agg(kkk,:)
2706 < cd                enddo
2707 < cd                do kkk=1,3
2708 < cd                write (iout,*) 'iii 2 kkk',kkk
2709 < cd                write (iout,*) aggi(kkk,:)
2710 < cd                enddo
2711 < cd                do kkk=1,3
2712 < cd                write (iout,*) 'iii 3 kkk',kkk
2713 < cd                write (iout,*) aggi1(kkk,:)
2714 < cd                enddo
2715 < cd                do kkk=1,3
2716 < cd                write (iout,*) 'iii 4 kkk',kkk
2717 < cd                write (iout,*) aggj(kkk,:)
2718 < cd                enddo
2719 < cd                do kkk=1,3
2720 < cd                write (iout,*) 'iii 5 kkk',kkk
2721 < cd                write (iout,*) aggj1(kkk,:)
2722 < cd                enddo
2723 2393,2395d3330
2724 < c                      do mm=1,5
2725 < c                      a_chuj_der(k,l,m,mm,num_conti,i)=0.0d0
2726 < c                      enddo
2727 2408,2409c3343,3356
2728 <                 ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
2729 <                 ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
2730 ---
2731 > c                 ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
2732 >                 ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
2733 >                 if (ees0tmp.gt.0) then
2734 >                   ees0pij=dsqrt(ees0tmp)
2735 >                 else
2736 >                   ees0pij=0
2737 >                 endif
2738 > c                ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
2739 >                 ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
2740 >                 if (ees0tmp.gt.0) then
2741 >                   ees0mij=dsqrt(ees0tmp)
2742 >                 else
2743 >                   ees0mij=0
2744 >                 endif
2745 2418,2421c3365,3366
2746 < c                write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij,
2747 < c     & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont
2748 <                 facont_hb(num_conti,i)=fcont
2749 <                 if (calc_grad) then
2750 ---
2751 > c               write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij,
2752 > c    & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont
2753 2448a3394
2754 >                 facont_hb(num_conti,i)=fcont
2755 2472,2474c3418,3424
2756 <                   ghalfp=0.5D0*gggp(k)
2757 <                   ghalfm=0.5D0*gggm(k)
2758 <                   gacontp_hb1(k,num_conti,i)=ghalfp
2759 ---
2760 > c
2761 > c 10/24/08 cgrad and ! comments indicate the parts of the code removed
2762 > c          following the change of gradient-summation algorithm.
2763 > c
2764 > cgrad                  ghalfp=0.5D0*gggp(k)
2765 > cgrad                  ghalfm=0.5D0*gggm(k)
2766 >                   gacontp_hb1(k,num_conti,i)=!ghalfp
2767 2477c3427
2768 <                   gacontp_hb2(k,num_conti,i)=ghalfp
2769 ---
2770 >                   gacontp_hb2(k,num_conti,i)=!ghalfp
2771 2481c3431
2772 <                   gacontm_hb1(k,num_conti,i)=ghalfm
2773 ---
2774 >                   gacontm_hb1(k,num_conti,i)=!ghalfm
2775 2484c3434
2776 <                   gacontm_hb2(k,num_conti,i)=ghalfm
2777 ---
2778 >                   gacontm_hb2(k,num_conti,i)=!ghalfm
2779 2489d3438
2780 <                 endif
2781 2503,2513c3452,3469
2782 <  1216     continue
2783 <         enddo ! j
2784 <         num_cont_hb(i)=num_conti
2785 <  1215   continue
2786 <       enddo   ! i
2787 < cd      do i=1,nres
2788 < cd        write (iout,'(i3,3f10.5,5x,3f10.5)')
2789 < cd     &     i,(gel_loc(k,i),k=1,3),gel_loc_loc(i)
2790 < cd      enddo
2791 < c 12/7/99 Adam eello_turn3 will be considered as a separate energy term
2792 < ccc      eel_loc=eel_loc+eello_turn3
2793 ---
2794 >           if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
2795 >             do k=1,4
2796 >               do l=1,3
2797 >                 ghalf=0.5d0*agg(l,k)
2798 >                 aggi(l,k)=aggi(l,k)+ghalf
2799 >                 aggi1(l,k)=aggi1(l,k)+agg(l,k)
2800 >                 aggj(l,k)=aggj(l,k)+ghalf
2801 >               enddo
2802 >             enddo
2803 >             if (j.eq.nres-1 .and. i.lt.j-2) then
2804 >               do k=1,4
2805 >                 do l=1,3
2806 >                   aggj1(l,k)=aggj1(l,k)+agg(l,k)
2807 >                 enddo
2808 >               enddo
2809 >             endif
2810 >           endif
2811 > c          t_eelecij=t_eelecij+MPI_Wtime()-time00
2812 2517c3473
2813 <       subroutine eturn34(i,j,eello_turn3,eello_turn4)
2814 ---
2815 >       subroutine eturn3(i,eello_turn3)
2816 2521d3476
2817 <       include 'DIMENSIONS.ZSCOPT'
2818 2532a3488
2819 >       include 'COMMON.CONTROL'
2820 2538,2540c3494,3503
2821 <      &    aggj(3,4),aggj1(3,4),a_temp(2,2)
2822 <       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,j1,j2
2823 <       if (j.eq.i+2) then
2824 ---
2825 >      &    aggj(3,4),aggj1(3,4),a_temp(2,2),auxmat3(2,2)
2826 >       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
2827 >      &    dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
2828 >      &    num_conti,j1,j2
2829 >       j=i+2
2830 > c      write (iout,*) "eturn3",i,j,j1,j2
2831 >       a_temp(1,1)=a22
2832 >       a_temp(1,2)=a23
2833 >       a_temp(2,1)=a32
2834 >       a_temp(2,2)=a33
2835 2555a3519,3520
2836 >         if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
2837 >      &          'eturn3',i,j,0.5d0*(pizda(1,1)+pizda(2,2))
2838 2559d3523
2839 <         if (calc_grad) then
2840 2562,2563c3526,3527
2841 <         call transpose2(auxmat2(1,1),pizda(1,1))
2842 <         call matmat2(a_temp(1,1),pizda(1,1),pizda(1,1))
2843 ---
2844 >         call transpose2(auxmat2(1,1),auxmat3(1,1))
2845 >         call matmat2(a_temp(1,1),auxmat3(1,1),pizda(1,1))
2846 2567,2568c3531,3532
2847 <         call transpose2(auxmat2(1,1),pizda(1,1))
2848 <         call matmat2(a_temp(1,1),pizda(1,1),pizda(1,1))
2849 ---
2850 >         call transpose2(auxmat2(1,1),auxmat3(1,1))
2851 >         call matmat2(a_temp(1,1),auxmat3(1,1),pizda(1,1))
2852 2573,2576c3537,3544
2853 <           a_temp(1,1)=aggi(l,1)
2854 <           a_temp(1,2)=aggi(l,2)
2855 <           a_temp(2,1)=aggi(l,3)
2856 <           a_temp(2,2)=aggi(l,4)
2857 ---
2858 > c            ghalf1=0.5d0*agg(l,1)
2859 > c            ghalf2=0.5d0*agg(l,2)
2860 > c            ghalf3=0.5d0*agg(l,3)
2861 > c            ghalf4=0.5d0*agg(l,4)
2862 >           a_temp(1,1)=aggi(l,1)!+ghalf1
2863 >           a_temp(1,2)=aggi(l,2)!+ghalf2
2864 >           a_temp(2,1)=aggi(l,3)!+ghalf3
2865 >           a_temp(2,2)=aggi(l,4)!+ghalf4
2866 2580,2583c3548,3551
2867 <           a_temp(1,1)=aggi1(l,1)
2868 <           a_temp(1,2)=aggi1(l,2)
2869 <           a_temp(2,1)=aggi1(l,3)
2870 <           a_temp(2,2)=aggi1(l,4)
2871 ---
2872 >           a_temp(1,1)=aggi1(l,1)!+agg(l,1)
2873 >           a_temp(1,2)=aggi1(l,2)!+agg(l,2)
2874 >           a_temp(2,1)=aggi1(l,3)!+agg(l,3)
2875 >           a_temp(2,2)=aggi1(l,4)!+agg(l,4)
2876 2587,2590c3555,3558
2877 <           a_temp(1,1)=aggj(l,1)
2878 <           a_temp(1,2)=aggj(l,2)
2879 <           a_temp(2,1)=aggj(l,3)
2880 <           a_temp(2,2)=aggj(l,4)
2881 ---
2882 >           a_temp(1,1)=aggj(l,1)!+ghalf1
2883 >           a_temp(1,2)=aggj(l,2)!+ghalf2
2884 >           a_temp(2,1)=aggj(l,3)!+ghalf3
2885 >           a_temp(2,2)=aggj(l,4)!+ghalf4
2886 2602,2603c3570,3598
2887 <         endif
2888 <       else if (j.eq.i+3 .and. itype(i+2).ne.21) then
2889 ---
2890 >       return
2891 >       end
2892 > C-------------------------------------------------------------------------------
2893 >       subroutine eturn4(i,eello_turn4)
2894 > C Third- and fourth-order contributions from turns
2895 >       implicit real*8 (a-h,o-z)
2896 >       include 'DIMENSIONS'
2897 >       include 'COMMON.IOUNITS'
2898 >       include 'COMMON.GEO'
2899 >       include 'COMMON.VAR'
2900 >       include 'COMMON.LOCAL'
2901 >       include 'COMMON.CHAIN'
2902 >       include 'COMMON.DERIV'
2903 >       include 'COMMON.INTERACT'
2904 >       include 'COMMON.CONTACTS'
2905 >       include 'COMMON.TORSION'
2906 >       include 'COMMON.VECTORS'
2907 >       include 'COMMON.FFIELD'
2908 >       include 'COMMON.CONTROL'
2909 >       dimension ggg(3)
2910 >       double precision auxmat(2,2),auxmat1(2,2),auxmat2(2,2),pizda(2,2),
2911 >      &  e1t(2,2),e2t(2,2),e3t(2,2),e1tder(2,2),e2tder(2,2),e3tder(2,2),
2912 >      &  e1a(2,2),ae3(2,2),ae3e2(2,2),auxvec(2),auxvec1(2)
2913 >       double precision agg(3,4),aggi(3,4),aggi1(3,4),
2914 >      &    aggj(3,4),aggj1(3,4),a_temp(2,2),auxmat3(2,2)
2915 >       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
2916 >      &    dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
2917 >      &    num_conti,j1,j2
2918 >       j=i+3
2919 2615a3611,3615
2920 > c        write (iout,*) "eturn4 i",i," j",j," j1",j1," j2",j2
2921 >         a_temp(1,1)=a22
2922 >         a_temp(1,2)=a23
2923 >         a_temp(2,1)=a32
2924 >         a_temp(2,2)=a33
2925 2618a3619
2926 > c        write(iout,*) "iti1",iti1," iti2",iti2," iti3",iti3
2927 2631a3633,3634
2928 >         if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
2929 >      &      'eturn4',i,j,-(s1+s2+s3)
2930 2635d3637
2931 <         if (calc_grad) then
2932 2658,2659c3660,3661
2933 <         call matmat2(auxmat(1,1),e2t(1,1),auxmat(1,1))
2934 <         call matmat2(auxmat(1,1),e1t(1,1),pizda(1,1))
2935 ---
2936 >         call matmat2(auxmat(1,1),e2t(1,1),auxmat3(1,1))
2937 >         call matmat2(auxmat3(1,1),e1t(1,1),pizda(1,1))
2938 2739a3742
2939 > c          write (iout,*) "s1",s1," s2",s2," s3",s3," s1+s2+s3",s1+s2+s3
2940 2742,2743d3744
2941 <         endif
2942 <       endif
2943 2779a3781,3877
2944 >       subroutine escp_soft_sphere(evdw2,evdw2_14)
2945 > C
2946 > C This subroutine calculates the excluded-volume interaction energy between
2947 > C peptide-group centers and side chains and its gradient in virtual-bond and
2948 > C side-chain vectors.
2949 > C
2950 >       implicit real*8 (a-h,o-z)
2951 >       include 'DIMENSIONS'
2952 >       include 'COMMON.GEO'
2953 >       include 'COMMON.VAR'
2954 >       include 'COMMON.LOCAL'
2955 >       include 'COMMON.CHAIN'
2956 >       include 'COMMON.DERIV'
2957 >       include 'COMMON.INTERACT'
2958 >       include 'COMMON.FFIELD'
2959 >       include 'COMMON.IOUNITS'
2960 >       include 'COMMON.CONTROL'
2961 >       dimension ggg(3)
2962 >       evdw2=0.0D0
2963 >       evdw2_14=0.0d0
2964 >       r0_scp=4.5d0
2965 > cd    print '(a)','Enter ESCP'
2966 > cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2967 >       do i=iatscp_s,iatscp_e
2968 >         if (itype(i).eq.21 .or. itype(i+1).eq.21) cycle
2969 >         iteli=itel(i)
2970 >         xi=0.5D0*(c(1,i)+c(1,i+1))
2971 >         yi=0.5D0*(c(2,i)+c(2,i+1))
2972 >         zi=0.5D0*(c(3,i)+c(3,i+1))
2973 >
2974 >         do iint=1,nscp_gr(i)
2975 >
2976 >         do j=iscpstart(i,iint),iscpend(i,iint)
2977 >           if (itype(j).eq.21) cycle
2978 >           itypj=itype(j)
2979 > C Uncomment following three lines for SC-p interactions
2980 > c         xj=c(1,nres+j)-xi
2981 > c         yj=c(2,nres+j)-yi
2982 > c         zj=c(3,nres+j)-zi
2983 > C Uncomment following three lines for Ca-p interactions
2984 >           xj=c(1,j)-xi
2985 >           yj=c(2,j)-yi
2986 >           zj=c(3,j)-zi
2987 >           rij=xj*xj+yj*yj+zj*zj
2988 >           r0ij=r0_scp
2989 >           r0ijsq=r0ij*r0ij
2990 >           if (rij.lt.r0ijsq) then
2991 >             evdwij=0.25d0*(rij-r0ijsq)**2
2992 >             fac=rij-r0ijsq
2993 >           else
2994 >             evdwij=0.0d0
2995 >             fac=0.0d0
2996 >           endif
2997 >           evdw2=evdw2+evdwij
2998 > C
2999 > C Calculate contributions to the gradient in the virtual-bond and SC vectors.
3000 > C
3001 >           ggg(1)=xj*fac
3002 >           ggg(2)=yj*fac
3003 >           ggg(3)=zj*fac
3004 > cgrad          if (j.lt.i) then
3005 > cd          write (iout,*) 'j<i'
3006 > C Uncomment following three lines for SC-p interactions
3007 > c           do k=1,3
3008 > c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
3009 > c           enddo
3010 > cgrad          else
3011 > cd          write (iout,*) 'j>i'
3012 > cgrad            do k=1,3
3013 > cgrad              ggg(k)=-ggg(k)
3014 > C Uncomment following line for SC-p interactions
3015 > c             gradx_scp(k,j)=gradx_scp(k,j)-ggg(k)
3016 > cgrad            enddo
3017 > cgrad          endif
3018 > cgrad          do k=1,3
3019 > cgrad            gvdwc_scp(k,i)=gvdwc_scp(k,i)-0.5D0*ggg(k)
3020 > cgrad          enddo
3021 > cgrad          kstart=min0(i+1,j)
3022 > cgrad          kend=max0(i-1,j-1)
3023 > cd        write (iout,*) 'i=',i,' j=',j,' kstart=',kstart,' kend=',kend
3024 > cd        write (iout,*) ggg(1),ggg(2),ggg(3)
3025 > cgrad          do k=kstart,kend
3026 > cgrad            do l=1,3
3027 > cgrad              gvdwc_scp(l,k)=gvdwc_scp(l,k)-ggg(l)
3028 > cgrad            enddo
3029 > cgrad          enddo
3030 >           do k=1,3
3031 >             gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
3032 >             gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
3033 >           enddo
3034 >         enddo
3035 >
3036 >         enddo ! iint
3037 >       enddo ! i
3038 >       return
3039 >       end
3040 > C-----------------------------------------------------------------------------
3041 2788d3885
3042 <       include 'DIMENSIONS.ZSCOPT'
3043 2796a3894
3044 >       include 'COMMON.CONTROL'
3045 2801,2802c3899
3046 < c      write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e,
3047 < c     &  ' scal14',scal14
3048 ---
3049 > cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
3050 2806,2808d3902
3051 < c        write (iout,*) "i",i," iteli",iteli," nscp_gr",nscp_gr(i),
3052 < c     &   " iscp",(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i))
3053 <         if (iteli.eq.0) goto 1225
3054 2836d3929
3055 < c          write (iout,*) i,j,evdwij
3056 2838c3931,3932
3057 <           if (calc_grad) then
3058 ---
3059 >           if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
3060 >      &        'evdw2',i,j,evdwij
3061 2846c3940
3062 <           if (j.lt.i) then
3063 ---
3064 > cgrad          if (j.lt.i) then
3065 2852c3946
3066 <           else
3067 ---
3068 > cgrad          else
3069 2854,2855c3948,3949
3070 <             do k=1,3
3071 <               ggg(k)=-ggg(k)
3072 ---
3073 > cgrad            do k=1,3
3074 > cgrad              ggg(k)=-ggg(k)
3075 2857,2864c3951,3959
3076 < c             gradx_scp(k,j)=gradx_scp(k,j)-ggg(k)
3077 <             enddo
3078 <           endif
3079 <           do k=1,3
3080 <             gvdwc_scp(k,i)=gvdwc_scp(k,i)-0.5D0*ggg(k)
3081 <           enddo
3082 <           kstart=min0(i+1,j)
3083 <           kend=max0(i-1,j-1)
3084 ---
3085 > ccgrad             gradx_scp(k,j)=gradx_scp(k,j)-ggg(k)
3086 > c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
3087 > cgrad            enddo
3088 > cgrad          endif
3089 > cgrad          do k=1,3
3090 > cgrad            gvdwc_scp(k,i)=gvdwc_scp(k,i)-0.5D0*ggg(k)
3091 > cgrad          enddo
3092 > cgrad          kstart=min0(i+1,j)
3093 > cgrad          kend=max0(i-1,j-1)
3094 2867,2870c3962,3969
3095 <           do k=kstart,kend
3096 <             do l=1,3
3097 <               gvdwc_scp(l,k)=gvdwc_scp(l,k)-ggg(l)
3098 <             enddo
3099 ---
3100 > cgrad          do k=kstart,kend
3101 > cgrad            do l=1,3
3102 > cgrad              gvdwc_scp(l,k)=gvdwc_scp(l,k)-ggg(l)
3103 > cgrad            enddo
3104 > cgrad          enddo
3105 >           do k=1,3
3106 >             gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
3107 >             gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
3108 2872d3970
3109 <           endif
3110 2873a3972
3111 >
3112 2875d3973
3113 <  1225   continue
3114 2879a3978
3115 >           gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
3116 2901d3999
3117 <       include 'DIMENSIONS.ZSCOPT'
3118 2906a4005
3119 >       include 'COMMON.IOUNITS'
3120 2909,2910c4008,4009
3121 < cd    print *,'edis: nhpb=',nhpb,' fbr=',fbr
3122 < cd    print *,'link_start=',link_start,' link_end=',link_end
3123 ---
3124 > cd      write(iout,*)'edis: nhpb=',nhpb,' fbr=',fbr
3125 > cd      write(iout,*)'link_start=',link_start,' link_end=',link_end
3126 2924a4024
3127 > cd        write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj
3128 2929a4030
3129 > cd          write (iout,*) "eij",eij
3130 2957,2960c4058,4065
3131 <         do j=iii,jjj-1
3132 <           do k=1,3
3133 <             ghpbc(k,j)=ghpbc(k,j)+ggg(k)
3134 <           enddo
3135 ---
3136 > cgrad        do j=iii,jjj-1
3137 > cgrad          do k=1,3
3138 > cgrad            ghpbc(k,j)=ghpbc(k,j)+ggg(k)
3139 > cgrad          enddo
3140 > cgrad        enddo
3141 >         do k=1,3
3142 >           ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
3143 >           ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
3144 2978d4082
3145 <       include 'DIMENSIONS.ZSCOPT'
3146 2994c4098,4099
3147 <       dsci_inv=dsc_inv(itypi)
3148 ---
3149 > c      dsci_inv=dsc_inv(itypi)
3150 >       dsci_inv=vbld_inv(nres+i)
3151 2996c4101,4102
3152 <       dscj_inv=dsc_inv(itypj)
3153 ---
3154 > c      dscj_inv=dsc_inv(itypj)
3155 >       dscj_inv=vbld_inv(nres+j)
3156 3034,3040c4140,4148
3157 <         gg(k)=ed*erij(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
3158 <       enddo
3159 <       do k=1,3
3160 <         ghpbx(k,i)=ghpbx(k,i)-gg(k)
3161 <      &            +(eom12*dc_norm(k,nres+j)+eom1*erij(k))*dsci_inv
3162 <         ghpbx(k,j)=ghpbx(k,j)+gg(k)
3163 <      &            +(eom12*dc_norm(k,nres+i)+eom2*erij(k))*dscj_inv
3164 ---
3165 >         ggk=ed*erij(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
3166 >         ghpbx(k,i)=ghpbx(k,i)-ggk
3167 >      &            +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
3168 >      &            +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
3169 >         ghpbx(k,j)=ghpbx(k,j)+ggk
3170 >      &            +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
3171 >      &            +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
3172 >         ghpbc(k,i)=ghpbc(k,i)-ggk
3173 >         ghpbc(k,j)=ghpbc(k,j)+ggk
3174 3045,3049c4153,4157
3175 <       do k=i,j-1
3176 <         do l=1,3
3177 <           ghpbc(l,k)=ghpbc(l,k)+gg(l)
3178 <         enddo
3179 <       enddo
3180 ---
3181 > cgrad      do k=i,j-1
3182 > cgrad        do l=1,3
3183 > cgrad          ghpbc(l,k)=ghpbc(l,k)+gg(l)
3184 > cgrad        enddo
3185 > cgrad      enddo
3186 3059d4166
3187 <       include 'DIMENSIONS.ZSCOPT'
3188 3070c4177
3189 <       logical energy_dec /.false./
3190 ---
3191 >       include 'COMMON.SETUP'
3192 3073,3074c4180,4181
3193 <       write (iout,*) "distchainmax",distchainmax
3194 <       do i=nnt+1,nct
3195 ---
3196 >       estr1=0.0d0
3197 >       do i=ibondp_start,ibondp_end
3198 3081c4188
3199 <           if (energy_dec) write(iout,*)
3200 ---
3201 >           if (energy_dec) write(iout,*)
3202 3085,3090c4192,4199
3203 <           diff = vbld(i)-vbldp0
3204 < c          write (iout,*) i,vbld(i),vbldp0,diff,AKP*diff*diff
3205 <           estr=estr+diff*diff
3206 <           do j=1,3
3207 <             gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i)
3208 <           enddo
3209 ---
3210 >         diff = vbld(i)-vbldp0
3211 >         if (energy_dec) write (iout,*)
3212 >      &     "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff
3213 >         estr=estr+diff*diff
3214 >         do j=1,3
3215 >           gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i)
3216 >         enddo
3217 > c        write (iout,'(i5,3f10.5)') i,(gradb(j,i-1),j=1,3)
3218 3092d4200
3219 <
3220 3094c4202
3221 <       estr=0.5d0*AKP*estr
3222 ---
3223 >       estr=0.5d0*AKP*estr+estr1
3224 3098c4206
3225 <       do i=nnt,nct
3226 ---
3227 >       do i=ibond_start,ibond_end
3228 3104,3105c4212,4214
3229 < c            write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff,
3230 < c     &      AKSC(1,iti),AKSC(1,iti)*diff*diff
3231 ---
3232 >             if (energy_dec) write (iout,*)
3233 >      &      "estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff,
3234 >      &      AKSC(1,iti),AKSC(1,iti)*diff*diff
3235 3112c4221
3236 <               diff=vbld(i+nres)-vbldsc0(j,iti)
3237 ---
3238 >               diff=vbld(i+nres)-vbldsc0(j,iti)
3239 3132c4241
3240 <               usumsqder=usumsqder+ud(j)*uprod2
3241 ---
3242 >               usumsqder=usumsqder+ud(j)*uprod2
3243 3134,3135d4242
3244 < c            write (iout,*) i,iti,vbld(i+nres),(vbldsc0(j,iti),
3245 < c     &      AKSC(j,iti),abond0(j,iti),u(j),j=1,nbi)
3246 3144c4251
3247 <       end
3248 ---
3249 >       end
3250 3154d4260
3251 <       include 'DIMENSIONS.ZSCOPT'
3252 3163a4270
3253 >       include 'COMMON.CONTROL'
3254 3169,3170c4276,4277
3255 <       time11=dexp(-2*time)
3256 <       time12=1.0d0
3257 ---
3258 > c      time11=dexp(-2*time)
3259 > c      time12=1.0d0
3260 3172d4278
3261 < c      write (iout,*) "nres",nres
3262 3174d4279
3263 < c      write (iout,*) ithet_start,ithet_end
3264 3182,3185c4287,4288
3265 <           phii=phi(i)
3266 <           icrc=0
3267 <           call proc_proc(phii,icrc)
3268 <           if (icrc.eq.1) phii=150.0
3269 ---
3270 >         phii=phi(i)
3271 >           if (phii.ne.phii) phii=150.0
3272 3191c4294
3273 <         else
3274 ---
3275 >         else
3276 3197,3200c4300,4301
3277 <           phii1=phi(i+1)
3278 <           icrc=0
3279 <           call proc_proc(phii1,icrc)
3280 <           if (icrc.eq.1) phii1=150.0
3281 ---
3282 >         phii1=phi(i+1)
3283 >           if (phii1.ne.phii1) phii1=150.0
3284 3211c4312
3285 <         endif
3286 ---
3287 >         endif
3288 3221d4321
3289 < c        write (iout,*) "thet_pred_mean",thet_pred_mean
3290 3224d4323
3291 < c        write (iout,*) "thet_pred_mean",thet_pred_mean
3292 3250,3251c4349,4350
3293 < c        write (iout,'(2i3,3f8.3,f10.5)') i,it,rad2deg*theta(i),
3294 < c     &    rad2deg*phii,rad2deg*phii1,ethetai
3295 ---
3296 >         if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
3297 >      &      'ebend',i,ethetai
3298 3255d4353
3299 <  1215   continue
3300 3379d4476
3301 <       include 'DIMENSIONS.ZSCOPT'
3302 3396d4492
3303 < c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
3304 3449,3451d4544
3305 < c        write (iout,*) "i",i," ityp1",itype(i-2),ityp1,
3306 < c     &   " ityp2",itype(i-1),ityp2," ityp3",itype(i),ityp3
3307 < c        call flush(iout)
3308 3571d4663
3309 <       include 'DIMENSIONS.ZSCOPT'
3310 3580a4673
3311 >       include 'COMMON.CONTROL'
3312 3598d4690
3313 < c        write (iout,*) "i",i," x",x(1),x(2),x(3)
3314 3672c4764,4766
3315 < c        write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc
3316 ---
3317 >         if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
3318 >      &     'escloc',i,escloci
3319 > c       write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc
3320 3745a4840,4844
3321 > #ifdef OSF
3322 >             adexp=bsc(j,it)-0.5D0*contr(j,iii)+emin
3323 >             if(adexp.ne.adexp) adexp=1.0
3324 >             expfac=dexp(adexp)
3325 > #else
3326 3746a4846
3327 > #endif
3328 3856d4955
3329 <       include 'DIMENSIONS.ZSCOPT'
3330 3991,3992c5090
3331 < c        write (2,*) "escloc",escloc
3332 <         if (.not. calc_grad) goto 1
3333 ---
3334 > c        write (2,*) "i",i," escloc",sumene,escloc
3335 4168a5267,5304
3336 > c------------------------------------------------------------------------------
3337 >       double precision function enesc(x,xx,yy,zz,cost2,sint2)
3338 >       implicit none
3339 >       double precision x(65),xx,yy,zz,cost2,sint2,sumene1,sumene2,
3340 >      & sumene3,sumene4,sumene,dsc_i,dp2_i,dscp1,dscp2,s1,s1_6,s2,s2_6
3341 >       sumene1= x(1)+  x(2)*xx+  x(3)*yy+  x(4)*zz+  x(5)*xx**2
3342 >      &   + x(6)*yy**2+  x(7)*zz**2+  x(8)*xx*zz+  x(9)*xx*yy
3343 >      &   + x(10)*yy*zz
3344 >       sumene2=  x(11) + x(12)*xx + x(13)*yy + x(14)*zz + x(15)*xx**2
3345 >      & + x(16)*yy**2 + x(17)*zz**2 + x(18)*xx*zz + x(19)*xx*yy
3346 >      & + x(20)*yy*zz
3347 >       sumene3=  x(21) +x(22)*xx +x(23)*yy +x(24)*zz +x(25)*xx**2
3348 >      &  +x(26)*yy**2 +x(27)*zz**2 +x(28)*xx*zz +x(29)*xx*yy
3349 >      &  +x(30)*yy*zz +x(31)*xx**3 +x(32)*yy**3 +x(33)*zz**3
3350 >      &  +x(34)*(xx**2)*yy +x(35)*(xx**2)*zz +x(36)*(yy**2)*xx
3351 >      &  +x(37)*(yy**2)*zz +x(38)*(zz**2)*xx +x(39)*(zz**2)*yy
3352 >      &  +x(40)*xx*yy*zz
3353 >       sumene4= x(41) +x(42)*xx +x(43)*yy +x(44)*zz +x(45)*xx**2
3354 >      &  +x(46)*yy**2 +x(47)*zz**2 +x(48)*xx*zz +x(49)*xx*yy
3355 >      &  +x(50)*yy*zz +x(51)*xx**3 +x(52)*yy**3 +x(53)*zz**3
3356 >      &  +x(54)*(xx**2)*yy +x(55)*(xx**2)*zz +x(56)*(yy**2)*xx
3357 >      &  +x(57)*(yy**2)*zz +x(58)*(zz**2)*xx +x(59)*(zz**2)*yy
3358 >      &  +x(60)*xx*yy*zz
3359 >       dsc_i   = 0.743d0+x(61)
3360 >       dp2_i   = 1.9d0+x(62)
3361 >       dscp1=dsqrt(dsc_i**2+dp2_i**2-2*dsc_i*dp2_i
3362 >      &          *(xx*cost2+yy*sint2))
3363 >       dscp2=dsqrt(dsc_i**2+dp2_i**2-2*dsc_i*dp2_i
3364 >      &          *(xx*cost2-yy*sint2))
3365 >       s1=(1+x(63))/(0.1d0 + dscp1)
3366 >       s1_6=(1+x(64))/(0.1d0 + dscp1**6)
3367 >       s2=(1+x(65))/(0.1d0 + dscp2)
3368 >       s2_6=(1+x(65))/(0.1d0 + dscp2**6)
3369 >       sumene = ( sumene3*sint2 + sumene1)*(s1+s1_6)
3370 >      & + (sumene4*cost2 +sumene2)*(s2+s2_6)
3371 >       enesc=sumene
3372 >       return
3373 >       end
3374 4207d5342
3375 <       include 'DIMENSIONS.ZSCOPT'
3376 4252c5387
3377 <       subroutine etor(etors,edihcnstr,fact)
3378 ---
3379 >       subroutine etor(etors,edihcnstr)
3380 4255d5389
3381 <       include 'DIMENSIONS.ZSCOPT'
3382 4266a5401
3383 >       include 'COMMON.CONTROL'
3384 4273c5408,5409
3385 <         if (itype(i-2).eq.21 .or. itype(i-1).eq.21
3386 ---
3387 >       etors_ii=0.0D0
3388 >         if (itype(i-2).eq.21 .or. itype(i-1).eq.21
3389 4286a5423
3390 >             if (energy_dec) etors_ii=etors_ii+etorsi-v1(1,3,3)
3391 4294a5432,5433
3392 >             if (energy_dec) etors_ii=etors_ii+
3393 >      &                              v2ij*sinphi+dabs(v1ij)+dabs(v2ij)
3394 4303a5443,5444
3395 >             if (energy_dec) etors_ii=etors_ii+
3396 >      &                  v1ij*cosphi+v2ij*sinphi+dabs(v1ij)+dabs(v2ij)
3397 4306a5448,5449
3398 >         if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
3399 >              'etor',i,etors_ii
3400 4311c5454
3401 <         gloc(i-3,icg)=gloc(i-3,icg)+wtor*fact*gloci
3402 ---
3403 >         gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci
3404 4335a5479,5483
3405 >       subroutine etor_d(etors_d)
3406 >       etors_d=0.0d0
3407 >       return
3408 >       end
3409 > c----------------------------------------------------------------------------
3410 4337c5485
3411 <       subroutine etor(etors,edihcnstr,fact)
3412 ---
3413 >       subroutine etor(etors,edihcnstr)
3414 4340d5487
3415 <       include 'DIMENSIONS.ZSCOPT'
3416 4351a5499
3417 >       include 'COMMON.CONTROL'
3418 4355c5503
3419 < c      lprn=.true.
3420 ---
3421 > c     lprn=.true.
3422 4358c5506
3423 <         if (itype(i-2).eq.21 .or. itype(i-1).eq.21
3424 ---
3425 >         if (itype(i-2).eq.21 .or. itype(i-1).eq.21
3426 4360c5508
3427 <         if (itel(i-2).eq.0 .or. itel(i-1).eq.0) goto 1215
3428 ---
3429 >         etors_ii=0.0D0
3430 4371a5520,5521
3431 >           if (energy_dec) etors_ii=etors_ii+
3432 >      &                v1ij*cosphi+v2ij*sinphi
3433 4387a5538,5539
3434 >           if (energy_dec) etors_ii=etors_ii+
3435 >      &                vl1ij*pom1
3436 4392a5545,5546
3437 >           if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
3438 >      &         'etor',i,etors_ii-v0(itori,itori1)
3439 4397c5551
3440 <         gloc(i-3,icg)=gloc(i-3,icg)+wtor*fact*gloci
3441 ---
3442 >         gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci
3443 4399d5552
3444 <  1215   continue
3445 4403c5556,5557
3446 <       do i=1,ndih_constr
3447 ---
3448 > c      do i=1,ndih_constr
3449 >       do i=idihconstr_start,idihconstr_end
3450 4407d5560
3451 <         edihi=0.0d0
3452 4412d5564
3453 <           edihi=0.25d0*ftors*difi**4
3454 4417d5568
3455 <           edihi=0.25d0*ftors*difi**4
3456 4419c5570
3457 <           difi=0.0d0
3458 ---
3459 >           difi=0.0
3460 4421,4424c5572,5574
3461 < c        write (iout,'(2i5,4f10.5,e15.5)') i,itori,phii,phi0(i),difi,
3462 < c     &    drange(i),edihi
3463 < !        write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii,
3464 < !     &    rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg)
3465 ---
3466 > cd        write (iout,'(2i5,4f8.3,2e14.5)') i,itori,rad2deg*phii,
3467 > cd     &    rad2deg*phi0(i),  rad2deg*drange(i),
3468 > cd     &    rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg)
3469 4426c5576
3470 < !      write (iout,*) 'edihcnstr',edihcnstr
3471 ---
3472 > cd       write (iout,*) 'edihcnstr',edihcnstr
3473 4430c5580
3474 <       subroutine etor_d(etors_d,fact2)
3475 ---
3476 >       subroutine etor_d(etors_d)
3477 4434d5583
3478 <       include 'DIMENSIONS.ZSCOPT'
3479 4451c5600
3480 <       do i=iphi_start,iphi_end-1
3481 ---
3482 >       do i=iphid_start,iphid_end
3483 4454,4455d5602
3484 <         if (itel(i-2).eq.0 .or. itel(i-1).eq.0 .or. itel(i).eq.0)
3485 <      &     goto 1215
3486 4496,4498c5643,5644
3487 <         gloc(i-3,icg)=gloc(i-3,icg)+wtor_d*fact2*gloci1
3488 <         gloc(i-2,icg)=gloc(i-2,icg)+wtor_d*fact2*gloci2
3489 <  1215   continue
3490 ---
3491 >         gloc(i-3,icg)=gloc(i-3,icg)+wtor_d*gloci1
3492 >         gloc(i-2,icg)=gloc(i-2,icg)+wtor_d*gloci2
3493 4509c5655
3494 < c        of residues computed from AM1 energy surfaces of terminally-blocked
3495 ---
3496 > c        of residues computed from AM1  energy surfaces of terminally-blocked
3497 4513d5658
3498 <       include 'DIMENSIONS.ZSCOPT'
3499 4551c5696
3500 <         gsccor_loc(i-3)=gloci
3501 ---
3502 >         gsccor_loc(i-3)=gsccor_loc(i-3)+gloci
3503 4555c5700
3504 < c------------------------------------------------------------------------------
3505 ---
3506 > c----------------------------------------------------------------------------
3507 4638,4702c5783,5794
3508 <       enddo
3509 <       do m=i,j-1
3510 <         do ll=1,3
3511 <           gradcorr(ll,m)=gradcorr(ll,m)+gx(ll)
3512 <         enddo
3513 <       enddo
3514 <       do m=k,l-1
3515 <         do ll=1,3
3516 <           gradcorr(ll,m)=gradcorr(ll,m)+gx1(ll)
3517 <         enddo
3518 <       enddo
3519 <       esccorr=-eij*ekl
3520 <       return
3521 <       end
3522 < c------------------------------------------------------------------------------
3523 < #ifdef MPL
3524 <       subroutine pack_buffer(dimen1,dimen2,atom,indx,buffer)
3525 <       implicit real*8 (a-h,o-z)
3526 <       include 'DIMENSIONS'
3527 <       integer dimen1,dimen2,atom,indx
3528 <       double precision buffer(dimen1,dimen2)
3529 <       double precision zapas
3530 <       common /contacts_hb/ zapas(3,20,maxres,7),
3531 <      &   facont_hb(20,maxres),ees0p(20,maxres),ees0m(20,maxres),
3532 <      &         num_cont_hb(maxres),jcont_hb(20,maxres)
3533 <       num_kont=num_cont_hb(atom)
3534 <       do i=1,num_kont
3535 <         do k=1,7
3536 <           do j=1,3
3537 <             buffer(i,indx+(k-1)*3+j)=zapas(j,i,atom,k)
3538 <           enddo ! j
3539 <         enddo ! k
3540 <         buffer(i,indx+22)=facont_hb(i,atom)
3541 <         buffer(i,indx+23)=ees0p(i,atom)
3542 <         buffer(i,indx+24)=ees0m(i,atom)
3543 <         buffer(i,indx+25)=dfloat(jcont_hb(i,atom))
3544 <       enddo ! i
3545 <       buffer(1,indx+26)=dfloat(num_kont)
3546 <       return
3547 <       end
3548 < c------------------------------------------------------------------------------
3549 <       subroutine unpack_buffer(dimen1,dimen2,atom,indx,buffer)
3550 <       implicit real*8 (a-h,o-z)
3551 <       include 'DIMENSIONS'
3552 <       integer dimen1,dimen2,atom,indx
3553 <       double precision buffer(dimen1,dimen2)
3554 <       double precision zapas
3555 <       common /contacts_hb/ zapas(3,20,maxres,7),
3556 <      &         facont_hb(20,maxres),ees0p(20,maxres),ees0m(20,maxres),
3557 <      &         num_cont_hb(maxres),jcont_hb(20,maxres)
3558 <       num_kont=buffer(1,indx+26)
3559 <       num_kont_old=num_cont_hb(atom)
3560 <       num_cont_hb(atom)=num_kont+num_kont_old
3561 <       do i=1,num_kont
3562 <         ii=i+num_kont_old
3563 <         do k=1,7
3564 <           do j=1,3
3565 <             zapas(j,ii,atom,k)=buffer(i,indx+(k-1)*3+j)
3566 <           enddo ! j
3567 <         enddo ! k
3568 <         facont_hb(ii,atom)=buffer(i,indx+22)
3569 <         ees0p(ii,atom)=buffer(i,indx+23)
3570 <         ees0m(ii,atom)=buffer(i,indx+24)
3571 <         jcont_hb(ii,atom)=buffer(i,indx+25)
3572 <       enddo ! i
3573 ---
3574 >       enddo
3575 >       do m=i,j-1
3576 >         do ll=1,3
3577 >           gradcorr(ll,m)=gradcorr(ll,m)+gx(ll)
3578 >         enddo
3579 >       enddo
3580 >       do m=k,l-1
3581 >         do ll=1,3
3582 >           gradcorr(ll,m)=gradcorr(ll,m)+gx1(ll)
3583 >         enddo
3584 >       enddo
3585 >       esccorr=-eij*ekl
3586 4706d5797
3587 < #endif
3588 4711d5801
3589 <       include 'DIMENSIONS.ZSCOPT'
3590 4713,4714c5803,5812
3591 < #ifdef MPL
3592 <       include 'COMMON.INFO'
3593 ---
3594 > #ifdef MPI
3595 >       include "mpif.h"
3596 >       parameter (max_cont=maxconts)
3597 >       parameter (max_dim=26)
3598 >       integer source,CorrelType,CorrelID,CorrelType1,CorrelID1,Error
3599 >       double precision zapas(max_dim,maxconts,max_fg_procs),
3600 >      &  zapas_recv(max_dim,maxconts,max_fg_procs)
3601 >       common /przechowalnia/ zapas
3602 >       integer status(MPI_STATUS_SIZE),req(maxconts*2),
3603 >      &  status_array(MPI_STATUS_SIZE,maxconts*2)
3604 4715a5814
3605 >       include 'COMMON.SETUP'
3606 4720,4728c5819,5821
3607 < #ifdef MPL
3608 <       parameter (max_cont=maxconts)
3609 <       parameter (max_dim=2*(8*3+2))
3610 <       parameter (msglen1=max_cont*max_dim*4)
3611 <       parameter (msglen2=2*msglen1)
3612 <       integer source,CorrelType,CorrelID,Error
3613 <       double precision buffer(max_cont,max_dim)
3614 < #endif
3615 <       double precision gx(3),gx1(3)
3616 ---
3617 >       include 'COMMON.CONTROL'
3618 >       include 'COMMON.LOCAL'
3619 >       double precision gx(3),gx1(3),time00
3620 4733c5826
3621 < #ifdef MPL
3622 ---
3623 > #ifdef MPI
3624 4736c5829
3625 <       if (fgProcs.le.1) goto 30
3626 ---
3627 >       if (nfgtasks.le.1) goto 30
3628 4738c5831
3629 <         write (iout,'(a)') 'Contact function values:'
3630 ---
3631 >         write (iout,'(a)') 'Contact function values before RECEIVE:'
3632 4745,4746c5838,5917
3633 < C Caution! Following code assumes that electrostatic interactions concerning
3634 < C a given atom are split among at most two processors!
3635 ---
3636 >       call flush(iout)
3637 >       do i=1,ntask_cont_from
3638 >         ncont_recv(i)=0
3639 >       enddo
3640 >       do i=1,ntask_cont_to
3641 >         ncont_sent(i)=0
3642 >       enddo
3643 > c      write (iout,*) "ntask_cont_from",ntask_cont_from," ntask_cont_to",
3644 > c     & ntask_cont_to
3645 > C Make the list of contacts to send to send to other procesors
3646 > c      write (iout,*) "limits",max0(iturn4_end-1,iatel_s),iturn3_end
3647 > c      call flush(iout)
3648 >       do i=iturn3_start,iturn3_end
3649 > c        write (iout,*) "make contact list turn3",i," num_cont",
3650 > c     &    num_cont_hb(i)
3651 >         call add_hb_contact(i,i+2,iturn3_sent_local(1,i))
3652 >       enddo
3653 >       do i=iturn4_start,iturn4_end
3654 > c        write (iout,*) "make contact list turn4",i," num_cont",
3655 > c     &   num_cont_hb(i)
3656 >         call add_hb_contact(i,i+3,iturn4_sent_local(1,i))
3657 >       enddo
3658 >       do ii=1,nat_sent
3659 >         i=iat_sent(ii)
3660 > c        write (iout,*) "make contact list longrange",i,ii," num_cont",
3661 > c     &    num_cont_hb(i)
3662 >         do j=1,num_cont_hb(i)
3663 >         do k=1,4
3664 >           jjc=jcont_hb(j,i)
3665 >           iproc=iint_sent_local(k,jjc,ii)
3666 > c          write (iout,*) "i",i," j",j," k",k," jjc",jjc," iproc",iproc
3667 >           if (iproc.gt.0) then
3668 >             ncont_sent(iproc)=ncont_sent(iproc)+1
3669 >             nn=ncont_sent(iproc)
3670 >             zapas(1,nn,iproc)=i
3671 >             zapas(2,nn,iproc)=jjc
3672 >             zapas(3,nn,iproc)=facont_hb(j,i)
3673 >             zapas(4,nn,iproc)=ees0p(j,i)
3674 >             zapas(5,nn,iproc)=ees0m(j,i)
3675 >             zapas(6,nn,iproc)=gacont_hbr(1,j,i)
3676 >             zapas(7,nn,iproc)=gacont_hbr(2,j,i)
3677 >             zapas(8,nn,iproc)=gacont_hbr(3,j,i)
3678 >             zapas(9,nn,iproc)=gacontm_hb1(1,j,i)
3679 >             zapas(10,nn,iproc)=gacontm_hb1(2,j,i)
3680 >             zapas(11,nn,iproc)=gacontm_hb1(3,j,i)
3681 >             zapas(12,nn,iproc)=gacontp_hb1(1,j,i)
3682 >             zapas(13,nn,iproc)=gacontp_hb1(2,j,i)
3683 >             zapas(14,nn,iproc)=gacontp_hb1(3,j,i)
3684 >             zapas(15,nn,iproc)=gacontm_hb2(1,j,i)
3685 >             zapas(16,nn,iproc)=gacontm_hb2(2,j,i)
3686 >             zapas(17,nn,iproc)=gacontm_hb2(3,j,i)
3687 >             zapas(18,nn,iproc)=gacontp_hb2(1,j,i)
3688 >             zapas(19,nn,iproc)=gacontp_hb2(2,j,i)
3689 >             zapas(20,nn,iproc)=gacontp_hb2(3,j,i)
3690 >             zapas(21,nn,iproc)=gacontm_hb3(1,j,i)
3691 >             zapas(22,nn,iproc)=gacontm_hb3(2,j,i)
3692 >             zapas(23,nn,iproc)=gacontm_hb3(3,j,i)
3693 >             zapas(24,nn,iproc)=gacontp_hb3(1,j,i)
3694 >             zapas(25,nn,iproc)=gacontp_hb3(2,j,i)
3695 >             zapas(26,nn,iproc)=gacontp_hb3(3,j,i)
3696 >           endif
3697 >         enddo
3698 >         enddo
3699 >       enddo
3700 >       if (lprn) then
3701 >       write (iout,*)
3702 >      &  "Numbers of contacts to be sent to other processors",
3703 >      &  (ncont_sent(i),i=1,ntask_cont_to)
3704 >       write (iout,*) "Contacts sent"
3705 >       do ii=1,ntask_cont_to
3706 >         nn=ncont_sent(ii)
3707 >         iproc=itask_cont_to(ii)
3708 >         write (iout,*) nn," contacts to processor",iproc,
3709 >      &   " of CONT_TO_COMM group"
3710 >         do i=1,nn
3711 >           write(iout,'(2f5.0,4f10.5)')(zapas(j,i,ii),j=1,5)
3712 >         enddo
3713 >       enddo
3714 >       call flush(iout)
3715 >       endif
3716 4748,4836c5919,6040
3717 <       CorrelID=MyID+1
3718 <       ldone=.false.
3719 <       do i=1,max_cont
3720 <         do j=1,max_dim
3721 <           buffer(i,j)=0.0D0
3722 <         enddo
3723 <       enddo
3724 <       mm=mod(MyRank,2)
3725 < cd    write (iout,*) 'MyRank',MyRank,' mm',mm
3726 <       if (mm) 20,20,10
3727 <    10 continue
3728 < cd    write (iout,*) 'Sending: MyRank',MyRank,' mm',mm,' ldone',ldone
3729 <       if (MyRank.gt.0) then
3730 < C Send correlation contributions to the preceding processor
3731 <         msglen=msglen1
3732 <         nn=num_cont_hb(iatel_s)
3733 <         call pack_buffer(max_cont,max_dim,iatel_s,0,buffer)
3734 < cd      write (iout,*) 'The BUFFER array:'
3735 < cd      do i=1,nn
3736 < cd        write (iout,'(i2,9(3f8.3,2x))') i,(buffer(i,j),j=1,26)
3737 < cd      enddo
3738 <         if (ielstart(iatel_s).gt.iatel_s+ispp) then
3739 <           msglen=msglen2
3740 <             call pack_buffer(max_cont,max_dim,iatel_s+1,26,buffer)
3741 < C Clear the contacts of the atom passed to the neighboring processor
3742 <         nn=num_cont_hb(iatel_s+1)
3743 < cd      do i=1,nn
3744 < cd        write (iout,'(i2,9(3f8.3,2x))') i,(buffer(i,j+26),j=1,26)
3745 < cd      enddo
3746 <             num_cont_hb(iatel_s)=0
3747 <         endif
3748 < cd      write (iout,*) 'Processor ',MyID,MyRank,
3749 < cd   & ' is sending correlation contribution to processor',MyID-1,
3750 < cd   & ' msglen=',msglen
3751 < cd      write (*,*) 'Processor ',MyID,MyRank,
3752 < cd   & ' is sending correlation contribution to processor',MyID-1,
3753 < cd   & ' msglen=',msglen,' CorrelType=',CorrelType
3754 <         call mp_bsend(buffer,msglen,MyID-1,CorrelType,CorrelID)
3755 < cd      write (iout,*) 'Processor ',MyID,
3756 < cd   & ' has sent correlation contribution to processor',MyID-1,
3757 < cd   & ' msglen=',msglen,' CorrelID=',CorrelID
3758 < cd      write (*,*) 'Processor ',MyID,
3759 < cd   & ' has sent correlation contribution to processor',MyID-1,
3760 < cd   & ' msglen=',msglen,' CorrelID=',CorrelID
3761 <         msglen=msglen1
3762 <       endif ! (MyRank.gt.0)
3763 <       if (ldone) goto 30
3764 <       ldone=.true.
3765 <    20 continue
3766 < cd    write (iout,*) 'Receiving: MyRank',MyRank,' mm',mm,' ldone',ldone
3767 <       if (MyRank.lt.fgProcs-1) then
3768 < C Receive correlation contributions from the next processor
3769 <         msglen=msglen1
3770 <         if (ielend(iatel_e).lt.nct-1) msglen=msglen2
3771 < cd      write (iout,*) 'Processor',MyID,
3772 < cd   & ' is receiving correlation contribution from processor',MyID+1,
3773 < cd   & ' msglen=',msglen,' CorrelType=',CorrelType
3774 < cd      write (*,*) 'Processor',MyID,
3775 < cd   & ' is receiving correlation contribution from processor',MyID+1,
3776 < cd   & ' msglen=',msglen,' CorrelType=',CorrelType
3777 <         nbytes=-1
3778 <         do while (nbytes.le.0)
3779 <           call mp_probe(MyID+1,CorrelType,nbytes)
3780 <         enddo
3781 < cd      print *,'Processor',MyID,' msglen',msglen,' nbytes',nbytes
3782 <         call mp_brecv(buffer,msglen,MyID+1,CorrelType,nbytes)
3783 < cd      write (iout,*) 'Processor',MyID,
3784 < cd   & ' has received correlation contribution from processor',MyID+1,
3785 < cd   & ' msglen=',msglen,' nbytes=',nbytes
3786 < cd      write (iout,*) 'The received BUFFER array:'
3787 < cd      do i=1,max_cont
3788 < cd        write (iout,'(i2,9(3f8.3,2x))') i,(buffer(i,j),j=1,52)
3789 < cd      enddo
3790 <         if (msglen.eq.msglen1) then
3791 <           call unpack_buffer(max_cont,max_dim,iatel_e+1,0,buffer)
3792 <         else if (msglen.eq.msglen2)  then
3793 <           call unpack_buffer(max_cont,max_dim,iatel_e,0,buffer)
3794 <           call unpack_buffer(max_cont,max_dim,iatel_e+1,26,buffer)
3795 <         else
3796 <           write (iout,*)
3797 <      & 'ERROR!!!! message length changed while processing correlations.'
3798 <           write (*,*)
3799 <      & 'ERROR!!!! message length changed while processing correlations.'
3800 <           call mp_stopall(Error)
3801 <         endif ! msglen.eq.msglen1
3802 <       endif ! MyRank.lt.fgProcs-1
3803 <       if (ldone) goto 30
3804 <       ldone=.true.
3805 <       goto 10
3806 ---
3807 >       CorrelID=fg_rank+1
3808 >       CorrelType1=478
3809 >       CorrelID1=nfgtasks+fg_rank+1
3810 >       ireq=0
3811 > C Receive the numbers of needed contacts from other processors
3812 >       do ii=1,ntask_cont_from
3813 >         iproc=itask_cont_from(ii)
3814 >         ireq=ireq+1
3815 >         call MPI_Irecv(ncont_recv(ii),1,MPI_INTEGER,iproc,CorrelType,
3816 >      &    FG_COMM,req(ireq),IERR)
3817 >       enddo
3818 > c      write (iout,*) "IRECV ended"
3819 > c      call flush(iout)
3820 > C Send the number of contacts needed by other processors
3821 >       do ii=1,ntask_cont_to
3822 >         iproc=itask_cont_to(ii)
3823 >         ireq=ireq+1
3824 >         call MPI_Isend(ncont_sent(ii),1,MPI_INTEGER,iproc,CorrelType,
3825 >      &    FG_COMM,req(ireq),IERR)
3826 >       enddo
3827 > c      write (iout,*) "ISEND ended"
3828 > c      write (iout,*) "number of requests (nn)",ireq
3829 >       call flush(iout)
3830 >       if (ireq.gt.0)
3831 >      &  call MPI_Waitall(ireq,req,status_array,ierr)
3832 > c      write (iout,*)
3833 > c     &  "Numbers of contacts to be received from other processors",
3834 > c     &  (ncont_recv(i),i=1,ntask_cont_from)
3835 > c      call flush(iout)
3836 > C Receive contacts
3837 >       ireq=0
3838 >       do ii=1,ntask_cont_from
3839 >         iproc=itask_cont_from(ii)
3840 >         nn=ncont_recv(ii)
3841 > c        write (iout,*) "Receiving",nn," contacts from processor",iproc,
3842 > c     &   " of CONT_TO_COMM group"
3843 >         call flush(iout)
3844 >         if (nn.gt.0) then
3845 >           ireq=ireq+1
3846 >           call MPI_Irecv(zapas_recv(1,1,ii),nn*max_dim,
3847 >      &    MPI_DOUBLE_PRECISION,iproc,CorrelType1,FG_COMM,req(ireq),IERR)
3848 > c          write (iout,*) "ireq,req",ireq,req(ireq)
3849 >         endif
3850 >       enddo
3851 > C Send the contacts to processors that need them
3852 >       do ii=1,ntask_cont_to
3853 >         iproc=itask_cont_to(ii)
3854 >         nn=ncont_sent(ii)
3855 > c        write (iout,*) nn," contacts to processor",iproc,
3856 > c     &   " of CONT_TO_COMM group"
3857 >         if (nn.gt.0) then
3858 >           ireq=ireq+1
3859 >           call MPI_Isend(zapas(1,1,ii),nn*max_dim,MPI_DOUBLE_PRECISION,
3860 >      &      iproc,CorrelType1,FG_COMM,req(ireq),IERR)
3861 > c          write (iout,*) "ireq,req",ireq,req(ireq)
3862 > c          do i=1,nn
3863 > c            write(iout,'(2f5.0,4f10.5)')(zapas(j,i,ii),j=1,5)
3864 > c          enddo
3865 >         endif
3866 >       enddo
3867 > c      write (iout,*) "number of requests (contacts)",ireq
3868 > c      write (iout,*) "req",(req(i),i=1,4)
3869 > c      call flush(iout)
3870 >       if (ireq.gt.0)
3871 >      & call MPI_Waitall(ireq,req,status_array,ierr)
3872 >       do iii=1,ntask_cont_from
3873 >         iproc=itask_cont_from(iii)
3874 >         nn=ncont_recv(iii)
3875 >         if (lprn) then
3876 >         write (iout,*) "Received",nn," contacts from processor",iproc,
3877 >      &   " of CONT_FROM_COMM group"
3878 >         call flush(iout)
3879 >         do i=1,nn
3880 >           write(iout,'(2f5.0,4f10.5)')(zapas_recv(j,i,iii),j=1,5)
3881 >         enddo
3882 >         call flush(iout)
3883 >         endif
3884 >         do i=1,nn
3885 >           ii=zapas_recv(1,i,iii)
3886 > c Flag the received contacts to prevent double-counting
3887 >           jj=-zapas_recv(2,i,iii)
3888 > c          write (iout,*) "iii",iii," i",i," ii",ii," jj",jj
3889 > c          call flush(iout)
3890 >           nnn=num_cont_hb(ii)+1
3891 >           num_cont_hb(ii)=nnn
3892 >           jcont_hb(nnn,ii)=jj
3893 >           facont_hb(nnn,ii)=zapas_recv(3,i,iii)
3894 >           ees0p(nnn,ii)=zapas_recv(4,i,iii)
3895 >           ees0m(nnn,ii)=zapas_recv(5,i,iii)
3896 >           gacont_hbr(1,nnn,ii)=zapas_recv(6,i,iii)
3897 >           gacont_hbr(2,nnn,ii)=zapas_recv(7,i,iii)
3898 >           gacont_hbr(3,nnn,ii)=zapas_recv(8,i,iii)
3899 >           gacontm_hb1(1,nnn,ii)=zapas_recv(9,i,iii)
3900 >           gacontm_hb1(2,nnn,ii)=zapas_recv(10,i,iii)
3901 >           gacontm_hb1(3,nnn,ii)=zapas_recv(11,i,iii)
3902 >           gacontp_hb1(1,nnn,ii)=zapas_recv(12,i,iii)
3903 >           gacontp_hb1(2,nnn,ii)=zapas_recv(13,i,iii)
3904 >           gacontp_hb1(3,nnn,ii)=zapas_recv(14,i,iii)
3905 >           gacontm_hb2(1,nnn,ii)=zapas_recv(15,i,iii)
3906 >           gacontm_hb2(2,nnn,ii)=zapas_recv(16,i,iii)
3907 >           gacontm_hb2(3,nnn,ii)=zapas_recv(17,i,iii)
3908 >           gacontp_hb2(1,nnn,ii)=zapas_recv(18,i,iii)
3909 >           gacontp_hb2(2,nnn,ii)=zapas_recv(19,i,iii)
3910 >           gacontp_hb2(3,nnn,ii)=zapas_recv(20,i,iii)
3911 >           gacontm_hb3(1,nnn,ii)=zapas_recv(21,i,iii)
3912 >           gacontm_hb3(2,nnn,ii)=zapas_recv(22,i,iii)
3913 >           gacontm_hb3(3,nnn,ii)=zapas_recv(23,i,iii)
3914 >           gacontp_hb3(1,nnn,ii)=zapas_recv(24,i,iii)
3915 >           gacontp_hb3(2,nnn,ii)=zapas_recv(25,i,iii)
3916 >           gacontp_hb3(3,nnn,ii)=zapas_recv(26,i,iii)
3917 >         enddo
3918 >       enddo
3919 >       call flush(iout)
3920 >       if (lprn) then
3921 >         write (iout,'(a)') 'Contact function values after receive:'
3922 >         do i=nnt,nct-2
3923 >           write (iout,'(2i3,50(1x,i3,f5.2))')
3924 >      &    i,num_cont_hb(i),(jcont_hb(j,i),facont_hb(j,i),
3925 >      &    j=1,num_cont_hb(i))
3926 >         enddo
3927 >         call flush(iout)
3928 >       endif
3929 4842c6046
3930 <           write (iout,'(2i3,50(1x,i2,f5.2))')
3931 ---
3932 >           write (iout,'(2i3,50(1x,i3,f5.2))')
3933 4856c6060
3934 <       do i=iatel_s,iatel_e+1
3935 ---
3936 >       do i=min0(iatel_s,iturn4_start),max0(iatel_e,iturn3_end)
3937 4861a6066
3938 >           jp=iabs(j)
3939 4863a6069
3940 >             jp1=iabs(j1)
3941 4866c6072,6074
3942 <             if (j1.eq.j+1 .or. j1.eq.j-1) then
3943 ---
3944 >             if ((j.gt.0 .and. j1.gt.0 .or. j.gt.0 .and. j1.lt.0
3945 >      &          .or. j.lt.0 .and. j1.gt.0) .and.
3946 >      &         (jp1.eq.jp+1 .or. jp1.eq.jp-1)) then
3947 4869c6077,6079
3948 <               ecorr=ecorr+ehbcorr(i,j,i+1,j1,jj,kk,0.72D0,0.32D0)
3949 ---
3950 >               ecorr=ecorr+ehbcorr(i,jp,i+1,jp1,jj,kk,0.72D0,0.32D0)
3951 >               if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
3952 >      &            'ecorrh',i,j,ehbcorr(i,j,i+1,j1,jj,kk,0.72D0,0.32D0)
3953 4891a6102,6158
3954 >       subroutine add_hb_contact(ii,jj,itask)
3955 >       implicit real*8 (a-h,o-z)
3956 >       include "DIMENSIONS"
3957 >       include "COMMON.IOUNITS"
3958 >       integer max_cont
3959 >       integer max_dim
3960 >       parameter (max_cont=maxconts)
3961 >       parameter (max_dim=26)
3962 >       include "COMMON.CONTACTS"
3963 >       double precision zapas(max_dim,maxconts,max_fg_procs),
3964 >      &  zapas_recv(max_dim,maxconts,max_fg_procs)
3965 >       common /przechowalnia/ zapas
3966 >       integer i,j,ii,jj,iproc,itask(4),nn
3967 > c      write (iout,*) "itask",itask
3968 >       do i=1,2
3969 >         iproc=itask(i)
3970 >         if (iproc.gt.0) then
3971 >           do j=1,num_cont_hb(ii)
3972 >             jjc=jcont_hb(j,ii)
3973 > c            write (iout,*) "i",ii," j",jj," jjc",jjc
3974 >             if (jjc.eq.jj) then
3975 >               ncont_sent(iproc)=ncont_sent(iproc)+1
3976 >               nn=ncont_sent(iproc)
3977 >               zapas(1,nn,iproc)=ii
3978 >               zapas(2,nn,iproc)=jjc
3979 >               zapas(3,nn,iproc)=facont_hb(j,ii)
3980 >               zapas(4,nn,iproc)=ees0p(j,ii)
3981 >               zapas(5,nn,iproc)=ees0m(j,ii)
3982 >               zapas(6,nn,iproc)=gacont_hbr(1,j,ii)
3983 >               zapas(7,nn,iproc)=gacont_hbr(2,j,ii)
3984 >               zapas(8,nn,iproc)=gacont_hbr(3,j,ii)
3985 >               zapas(9,nn,iproc)=gacontm_hb1(1,j,ii)
3986 >               zapas(10,nn,iproc)=gacontm_hb1(2,j,ii)
3987 >               zapas(11,nn,iproc)=gacontm_hb1(3,j,ii)
3988 >               zapas(12,nn,iproc)=gacontp_hb1(1,j,ii)
3989 >               zapas(13,nn,iproc)=gacontp_hb1(2,j,ii)
3990 >               zapas(14,nn,iproc)=gacontp_hb1(3,j,ii)
3991 >               zapas(15,nn,iproc)=gacontm_hb2(1,j,ii)
3992 >               zapas(16,nn,iproc)=gacontm_hb2(2,j,ii)
3993 >               zapas(17,nn,iproc)=gacontm_hb2(3,j,ii)
3994 >               zapas(18,nn,iproc)=gacontp_hb2(1,j,ii)
3995 >               zapas(19,nn,iproc)=gacontp_hb2(2,j,ii)
3996 >               zapas(20,nn,iproc)=gacontp_hb2(3,j,ii)
3997 >               zapas(21,nn,iproc)=gacontm_hb3(1,j,ii)
3998 >               zapas(22,nn,iproc)=gacontm_hb3(2,j,ii)
3999 >               zapas(23,nn,iproc)=gacontm_hb3(3,j,ii)
4000 >               zapas(24,nn,iproc)=gacontp_hb3(1,j,ii)
4001 >               zapas(25,nn,iproc)=gacontp_hb3(2,j,ii)
4002 >               zapas(26,nn,iproc)=gacontp_hb3(3,j,ii)
4003 >               exit
4004 >             endif
4005 >           enddo
4006 >         endif
4007 >       enddo
4008 >       return
4009 >       end
4010 > c------------------------------------------------------------------------------
4011 4897d6163
4012 <       include 'DIMENSIONS.ZSCOPT'
4013 4899,4900c6165,6174
4014 < #ifdef MPL
4015 <       include 'COMMON.INFO'
4016 ---
4017 > #ifdef MPI
4018 >       include "mpif.h"
4019 >       parameter (max_cont=maxconts)
4020 >       parameter (max_dim=70)
4021 >       integer source,CorrelType,CorrelID,CorrelType1,CorrelID1,Error
4022 >       double precision zapas(max_dim,maxconts,max_fg_procs),
4023 >      &  zapas_recv(max_dim,maxconts,max_fg_procs)
4024 >       common /przechowalnia/ zapas
4025 >       integer status(MPI_STATUS_SIZE),req(maxconts*2),
4026 >      &  status_array(MPI_STATUS_SIZE,maxconts*2)
4027 4901a6176
4028 >       include 'COMMON.SETUP'
4029 4903a6179
4030 >       include 'COMMON.LOCAL'
4031 4906,4913c6182,6183
4032 < #ifdef MPL
4033 <       parameter (max_cont=maxconts)
4034 <       parameter (max_dim=2*(8*3+2))
4035 <       parameter (msglen1=max_cont*max_dim*4)
4036 <       parameter (msglen2=2*msglen1)
4037 <       integer source,CorrelType,CorrelID,Error
4038 <       double precision buffer(max_cont,max_dim)
4039 < #endif
4040 ---
4041 >       include 'COMMON.CHAIN'
4042 >       include 'COMMON.CONTROL'
4043 4914a6185
4044 >       integer num_cont_hb_old(maxres)
4045 4916c6187,6188
4046 <
4047 ---
4048 >       double precision eello4,eello5,eelo6,eello_turn6
4049 >       external eello4,eello5,eello6,eello_turn6
4050 4920c6192,6195
4051 < #ifdef MPL
4052 ---
4053 > #ifdef MPI
4054 >       do i=1,nres
4055 >         num_cont_hb_old(i)=num_cont_hb(i)
4056 >       enddo
4057 4923c6198
4058 <       if (fgProcs.le.1) goto 30
4059 ---
4060 >       if (nfgtasks.le.1) goto 30
4061 4925c6200
4062 <         write (iout,'(a)') 'Contact function values:'
4063 ---
4064 >         write (iout,'(a)') 'Contact function values before RECEIVE:'
4065 4932,4933c6207,6282
4066 < C Caution! Following code assumes that electrostatic interactions concerning
4067 < C a given atom are split among at most two processors!
4068 ---
4069 >       call flush(iout)
4070 >       do i=1,ntask_cont_from
4071 >         ncont_recv(i)=0
4072 >       enddo
4073 >       do i=1,ntask_cont_to
4074 >         ncont_sent(i)=0
4075 >       enddo
4076 > c      write (iout,*) "ntask_cont_from",ntask_cont_from," ntask_cont_to",
4077 > c     & ntask_cont_to
4078 > C Make the list of contacts to send to send to other procesors
4079 >       do i=iturn3_start,iturn3_end
4080 > c        write (iout,*) "make contact list turn3",i," num_cont",
4081 > c     &    num_cont_hb(i)
4082 >         call add_hb_contact_eello(i,i+2,iturn3_sent_local(1,i))
4083 >       enddo
4084 >       do i=iturn4_start,iturn4_end
4085 > c        write (iout,*) "make contact list turn4",i," num_cont",
4086 > c     &   num_cont_hb(i)
4087 >         call add_hb_contact_eello(i,i+3,iturn4_sent_local(1,i))
4088 >       enddo
4089 >       do ii=1,nat_sent
4090 >         i=iat_sent(ii)
4091 > c        write (iout,*) "make contact list longrange",i,ii," num_cont",
4092 > c     &    num_cont_hb(i)
4093 >         do j=1,num_cont_hb(i)
4094 >         do k=1,4
4095 >           jjc=jcont_hb(j,i)
4096 >           iproc=iint_sent_local(k,jjc,ii)
4097 > c          write (iout,*) "i",i," j",j," k",k," jjc",jjc," iproc",iproc
4098 >           if (iproc.ne.0) then
4099 >             ncont_sent(iproc)=ncont_sent(iproc)+1
4100 >             nn=ncont_sent(iproc)
4101 >             zapas(1,nn,iproc)=i
4102 >             zapas(2,nn,iproc)=jjc
4103 >             zapas(3,nn,iproc)=d_cont(j,i)
4104 >             ind=3
4105 >             do kk=1,3
4106 >               ind=ind+1
4107 >               zapas(ind,nn,iproc)=grij_hb_cont(kk,j,i)
4108 >             enddo
4109 >             do kk=1,2
4110 >               do ll=1,2
4111 >                 ind=ind+1
4112 >                 zapas(ind,nn,iproc)=a_chuj(ll,kk,j,i)
4113 >               enddo
4114 >             enddo
4115 >             do jj=1,5
4116 >               do kk=1,3
4117 >                 do ll=1,2
4118 >                   do mm=1,2
4119 >                     ind=ind+1
4120 >                     zapas(ind,nn,iproc)=a_chuj_der(mm,ll,kk,jj,j,i)
4121 >                   enddo
4122 >                 enddo
4123 >               enddo
4124 >             enddo
4125 >           endif
4126 >         enddo
4127 >         enddo
4128 >       enddo
4129 >       if (lprn) then
4130 >       write (iout,*)
4131 >      &  "Numbers of contacts to be sent to other processors",
4132 >      &  (ncont_sent(i),i=1,ntask_cont_to)
4133 >       write (iout,*) "Contacts sent"
4134 >       do ii=1,ntask_cont_to
4135 >         nn=ncont_sent(ii)
4136 >         iproc=itask_cont_to(ii)
4137 >         write (iout,*) nn," contacts to processor",iproc,
4138 >      &   " of CONT_TO_COMM group"
4139 >         do i=1,nn
4140 >           write(iout,'(2f5.0,10f10.5)')(zapas(j,i,ii),j=1,10)
4141 >         enddo
4142 >       enddo
4143 >       call flush(iout)
4144 >       endif
4145 4935,5023c6284,6403
4146 <       CorrelID=MyID+1
4147 <       ldone=.false.
4148 <       do i=1,max_cont
4149 <         do j=1,max_dim
4150 <           buffer(i,j)=0.0D0
4151 <         enddo
4152 <       enddo
4153 <       mm=mod(MyRank,2)
4154 < cd    write (iout,*) 'MyRank',MyRank,' mm',mm
4155 <       if (mm) 20,20,10
4156 <    10 continue
4157 < cd    write (iout,*) 'Sending: MyRank',MyRank,' mm',mm,' ldone',ldone
4158 <       if (MyRank.gt.0) then
4159 < C Send correlation contributions to the preceding processor
4160 <         msglen=msglen1
4161 <         nn=num_cont_hb(iatel_s)
4162 <         call pack_buffer(max_cont,max_dim,iatel_s,0,buffer)
4163 < cd      write (iout,*) 'The BUFFER array:'
4164 < cd      do i=1,nn
4165 < cd        write (iout,'(i2,9(3f8.3,2x))') i,(buffer(i,j),j=1,26)
4166 < cd      enddo
4167 <         if (ielstart(iatel_s).gt.iatel_s+ispp) then
4168 <           msglen=msglen2
4169 <             call pack_buffer(max_cont,max_dim,iatel_s+1,26,buffer)
4170 < C Clear the contacts of the atom passed to the neighboring processor
4171 <         nn=num_cont_hb(iatel_s+1)
4172 < cd      do i=1,nn
4173 < cd        write (iout,'(i2,9(3f8.3,2x))') i,(buffer(i,j+26),j=1,26)
4174 < cd      enddo
4175 <             num_cont_hb(iatel_s)=0
4176 <         endif
4177 < cd      write (iout,*) 'Processor ',MyID,MyRank,
4178 < cd   & ' is sending correlation contribution to processor',MyID-1,
4179 < cd   & ' msglen=',msglen
4180 < cd      write (*,*) 'Processor ',MyID,MyRank,
4181 < cd   & ' is sending correlation contribution to processor',MyID-1,
4182 < cd   & ' msglen=',msglen,' CorrelType=',CorrelType
4183 <         call mp_bsend(buffer,msglen,MyID-1,CorrelType,CorrelID)
4184 < cd      write (iout,*) 'Processor ',MyID,
4185 < cd   & ' has sent correlation contribution to processor',MyID-1,
4186 < cd   & ' msglen=',msglen,' CorrelID=',CorrelID
4187 < cd      write (*,*) 'Processor ',MyID,
4188 < cd   & ' has sent correlation contribution to processor',MyID-1,
4189 < cd   & ' msglen=',msglen,' CorrelID=',CorrelID
4190 <         msglen=msglen1
4191 <       endif ! (MyRank.gt.0)
4192 <       if (ldone) goto 30
4193 <       ldone=.true.
4194 <    20 continue
4195 < cd    write (iout,*) 'Receiving: MyRank',MyRank,' mm',mm,' ldone',ldone
4196 <       if (MyRank.lt.fgProcs-1) then
4197 < C Receive correlation contributions from the next processor
4198 <         msglen=msglen1
4199 <         if (ielend(iatel_e).lt.nct-1) msglen=msglen2
4200 < cd      write (iout,*) 'Processor',MyID,
4201 < cd   & ' is receiving correlation contribution from processor',MyID+1,
4202 < cd   & ' msglen=',msglen,' CorrelType=',CorrelType
4203 < cd      write (*,*) 'Processor',MyID,
4204 < cd   & ' is receiving correlation contribution from processor',MyID+1,
4205 < cd   & ' msglen=',msglen,' CorrelType=',CorrelType
4206 <         nbytes=-1
4207 <         do while (nbytes.le.0)
4208 <           call mp_probe(MyID+1,CorrelType,nbytes)
4209 <         enddo
4210 < cd      print *,'Processor',MyID,' msglen',msglen,' nbytes',nbytes
4211 <         call mp_brecv(buffer,msglen,MyID+1,CorrelType,nbytes)
4212 < cd      write (iout,*) 'Processor',MyID,
4213 < cd   & ' has received correlation contribution from processor',MyID+1,
4214 < cd   & ' msglen=',msglen,' nbytes=',nbytes
4215 < cd      write (iout,*) 'The received BUFFER array:'
4216 < cd      do i=1,max_cont
4217 < cd        write (iout,'(i2,9(3f8.3,2x))') i,(buffer(i,j),j=1,52)
4218 < cd      enddo
4219 <         if (msglen.eq.msglen1) then
4220 <           call unpack_buffer(max_cont,max_dim,iatel_e+1,0,buffer)
4221 <         else if (msglen.eq.msglen2)  then
4222 <           call unpack_buffer(max_cont,max_dim,iatel_e,0,buffer)
4223 <           call unpack_buffer(max_cont,max_dim,iatel_e+1,26,buffer)
4224 <         else
4225 <           write (iout,*)
4226 <      & 'ERROR!!!! message length changed while processing correlations.'
4227 <           write (*,*)
4228 <      & 'ERROR!!!! message length changed while processing correlations.'
4229 <           call mp_stopall(Error)
4230 <         endif ! msglen.eq.msglen1
4231 <       endif ! MyRank.lt.fgProcs-1
4232 <       if (ldone) goto 30
4233 <       ldone=.true.
4234 <       goto 10
4235 ---
4236 >       CorrelID=fg_rank+1
4237 >       CorrelType1=478
4238 >       CorrelID1=nfgtasks+fg_rank+1
4239 >       ireq=0
4240 > C Receive the numbers of needed contacts from other processors
4241 >       do ii=1,ntask_cont_from
4242 >         iproc=itask_cont_from(ii)
4243 >         ireq=ireq+1
4244 >         call MPI_Irecv(ncont_recv(ii),1,MPI_INTEGER,iproc,CorrelType,
4245 >      &    FG_COMM,req(ireq),IERR)
4246 >       enddo
4247 > c      write (iout,*) "IRECV ended"
4248 > c      call flush(iout)
4249 > C Send the number of contacts needed by other processors
4250 >       do ii=1,ntask_cont_to
4251 >         iproc=itask_cont_to(ii)
4252 >         ireq=ireq+1
4253 >         call MPI_Isend(ncont_sent(ii),1,MPI_INTEGER,iproc,CorrelType,
4254 >      &    FG_COMM,req(ireq),IERR)
4255 >       enddo
4256 > c      write (iout,*) "ISEND ended"
4257 > c      write (iout,*) "number of requests (nn)",ireq
4258 >       call flush(iout)
4259 >       if (ireq.gt.0)
4260 >      &  call MPI_Waitall(ireq,req,status_array,ierr)
4261 > c      write (iout,*)
4262 > c     &  "Numbers of contacts to be received from other processors",
4263 > c     &  (ncont_recv(i),i=1,ntask_cont_from)
4264 > c      call flush(iout)
4265 > C Receive contacts
4266 >       ireq=0
4267 >       do ii=1,ntask_cont_from
4268 >         iproc=itask_cont_from(ii)
4269 >         nn=ncont_recv(ii)
4270 > c        write (iout,*) "Receiving",nn," contacts from processor",iproc,
4271 > c     &   " of CONT_TO_COMM group"
4272 >         call flush(iout)
4273 >         if (nn.gt.0) then
4274 >           ireq=ireq+1
4275 >           call MPI_Irecv(zapas_recv(1,1,ii),nn*max_dim,
4276 >      &    MPI_DOUBLE_PRECISION,iproc,CorrelType1,FG_COMM,req(ireq),IERR)
4277 > c          write (iout,*) "ireq,req",ireq,req(ireq)
4278 >         endif
4279 >       enddo
4280 > C Send the contacts to processors that need them
4281 >       do ii=1,ntask_cont_to
4282 >         iproc=itask_cont_to(ii)
4283 >         nn=ncont_sent(ii)
4284 > c        write (iout,*) nn," contacts to processor",iproc,
4285 > c     &   " of CONT_TO_COMM group"
4286 >         if (nn.gt.0) then
4287 >           ireq=ireq+1
4288 >           call MPI_Isend(zapas(1,1,ii),nn*max_dim,MPI_DOUBLE_PRECISION,
4289 >      &      iproc,CorrelType1,FG_COMM,req(ireq),IERR)
4290 > c          write (iout,*) "ireq,req",ireq,req(ireq)
4291 > c          do i=1,nn
4292 > c            write(iout,'(2f5.0,4f10.5)')(zapas(j,i,ii),j=1,5)
4293 > c          enddo
4294 >         endif
4295 >       enddo
4296 > c      write (iout,*) "number of requests (contacts)",ireq
4297 > c      write (iout,*) "req",(req(i),i=1,4)
4298 > c      call flush(iout)
4299 >       if (ireq.gt.0)
4300 >      & call MPI_Waitall(ireq,req,status_array,ierr)
4301 >       do iii=1,ntask_cont_from
4302 >         iproc=itask_cont_from(iii)
4303 >         nn=ncont_recv(iii)
4304 >         if (lprn) then
4305 >         write (iout,*) "Received",nn," contacts from processor",iproc,
4306 >      &   " of CONT_FROM_COMM group"
4307 >         call flush(iout)
4308 >         do i=1,nn
4309 >           write(iout,'(2f5.0,10f10.5)')(zapas_recv(j,i,iii),j=1,10)
4310 >         enddo
4311 >         call flush(iout)
4312 >         endif
4313 >         do i=1,nn
4314 >           ii=zapas_recv(1,i,iii)
4315 > c Flag the received contacts to prevent double-counting
4316 >           jj=-zapas_recv(2,i,iii)
4317 > c          write (iout,*) "iii",iii," i",i," ii",ii," jj",jj
4318 > c          call flush(iout)
4319 >           nnn=num_cont_hb(ii)+1
4320 >           num_cont_hb(ii)=nnn
4321 >           jcont_hb(nnn,ii)=jj
4322 >           d_cont(nnn,ii)=zapas_recv(3,i,iii)
4323 >           ind=3
4324 >           do kk=1,3
4325 >             ind=ind+1
4326 >             grij_hb_cont(kk,nnn,ii)=zapas_recv(ind,i,iii)
4327 >           enddo
4328 >           do kk=1,2
4329 >             do ll=1,2
4330 >               ind=ind+1
4331 >               a_chuj(ll,kk,nnn,ii)=zapas_recv(ind,i,iii)
4332 >             enddo
4333 >           enddo
4334 >           do jj=1,5
4335 >             do kk=1,3
4336 >               do ll=1,2
4337 >                 do mm=1,2
4338 >                   ind=ind+1
4339 >                   a_chuj_der(mm,ll,kk,jj,nnn,ii)=zapas_recv(ind,i,iii)
4340 >                 enddo
4341 >               enddo
4342 >             enddo
4343 >           enddo
4344 >         enddo
4345 >       enddo
4346 >       call flush(iout)
4347 >       if (lprn) then
4348 >         write (iout,'(a)') 'Contact function values after receive:'
4349 >         do i=nnt,nct-2
4350 >           write (iout,'(2i3,50(1x,i3,5f6.3))')
4351 >      &    i,num_cont_hb(i),(jcont_hb(j,i),d_cont(j,i),
4352 >      &    ((a_chuj(ll,kk,j,i),ll=1,2),kk=1,2),j=1,num_cont_hb(i))
4353 >         enddo
4354 >         call flush(iout)
4355 >       endif
4356 5029,5031c6409,6411
4357 <           write (iout,'(2i3,50(1x,i2,f5.2))')
4358 <      &    i,num_cont_hb(i),(jcont_hb(j,i),facont_hb(j,i),
4359 <      &    j=1,num_cont_hb(i))
4360 ---
4361 >           write (iout,'(2i3,50(1x,i2,5f6.3))')
4362 >      &    i,num_cont_hb(i),(jcont_hb(j,i),d_cont(j,i),
4363 >      &    ((a_chuj(ll,kk,j,i),ll=1,2),kk=1,2),j=1,num_cont_hb(i))
4364 5049a6430
4365 > #ifdef MOMENT
4366 5050a6432
4367 > #endif
4368 5055c6437,6443
4369 <       do i=iatel_s,iatel_e+1
4370 ---
4371 > c                write (iout,*) "gradcorr5 in eello5 before loop"
4372 > c                do iii=1,nres
4373 > c                  write (iout,'(i5,3f10.5)')
4374 > c     &             iii,(gradcorr5(jjj,iii),jjj=1,3)
4375 > c                enddo
4376 >       do i=min0(iatel_s,iturn4_start),max0(iatel_e+1,iturn3_end+1)
4377 > c        write (iout,*) "corr loop i",i
4378 5060a6449
4379 >           jp=iabs(j)
4380 5063c6452,6453
4381 < c            write (*,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1,
4382 ---
4383 >             jp1=iabs(j1)
4384 > c            write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1,
4385 5065c6455,6458
4386 <             if (j1.eq.j+1 .or. j1.eq.j-1) then
4387 ---
4388 > c            if (j1.eq.j+1 .or. j1.eq.j-1) then
4389 >             if ((j.gt.0 .and. j1.gt.0 .or. j.gt.0 .and. j1.lt.0
4390 >      &          .or. j.lt.0 .and. j1.gt.0) .and.
4391 >      &         (jp1.eq.jp+1 .or. jp1.eq.jp-1)) then
4392 5075,5076c6468,6469
4393 < c               write (*,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1,
4394 < c     &         ' jj=',jj,' kk=',kk
4395 ---
4396 > cd               write (iout,*) 'i=',i,' j=',jp,' i1=',i1,' j1=',jp1,
4397 > cd     &         ' jj=',jj,' kk=',kk
4398 5085,5086c6478,6483
4399 < cd     &          ' ekont=',ekont,' fprim=',fprimcont
4400 <                 call calc_eello(i,j,i+1,j1,jj,kk)
4401 ---
4402 > cd     &          ' ekont=',ekont,' fprim=',fprimcont,
4403 > cd     &          ' fac_prim1',fac_prim1,' fac_prim2',fac_prim2
4404 > cd               write (iout,*) "g_contij",g_contij
4405 > cd               write (iout,*) "grij_hb_cont i",grij_hb_cont(:,jj,i)
4406 > cd               write (iout,*) "grij_hb_cont i1",grij_hb_cont(:,jj,i1)
4407 >                 call calc_eello(i,jp,i+1,jp1,jj,kk)
4408 5088c6485,6493
4409 <      &            ecorr=ecorr+eello4(i,j,i+1,j1,jj,kk)
4410 ---
4411 >      &            ecorr=ecorr+eello4(i,jp,i+1,jp1,jj,kk)
4412 >                   if (energy_dec.and.wcorr4.gt.0.0d0)
4413 >      1                 write (iout,'(a6,4i5,0pf7.3)')
4414 >      2                'ecorr4',i,j,i+1,j1,eello4(i,jp,i+1,jp1,jj,kk)
4415 > c                write (iout,*) "gradcorr5 before eello5"
4416 > c                do iii=1,nres
4417 > c                  write (iout,'(i5,3f10.5)')
4418 > c     &             iii,(gradcorr5(jjj,iii),jjj=1,3)
4419 > c                enddo
4420 5090,5091c6495,6503
4421 <      &            ecorr5=ecorr5+eello5(i,j,i+1,j1,jj,kk)
4422 < c                print *,"wcorr5",ecorr5
4423 ---
4424 >      &            ecorr5=ecorr5+eello5(i,jp,i+1,jp1,jj,kk)
4425 > c                write (iout,*) "gradcorr5 after eello5"
4426 > c                do iii=1,nres
4427 > c                  write (iout,'(i5,3f10.5)')
4428 > c     &             iii,(gradcorr5(jjj,iii),jjj=1,3)
4429 > c                enddo
4430 >                   if (energy_dec.and.wcorr5.gt.0.0d0)
4431 >      1                 write (iout,'(a6,4i5,0pf7.3)')
4432 >      2                'ecorr5',i,j,i+1,j1,eello5(i,jp,i+1,jp1,jj,kk)
4433 5093,5094c6505,6506
4434 < cd                write(2,*)'ijkl',i,j,i+1,j1
4435 <                 if (wcorr6.gt.0.0d0 .and. (j.ne.i+4 .or. j1.ne.i+3
4436 ---
4437 > cd                write(2,*)'ijkl',i,jp,i+1,jp1
4438 >                 if (wcorr6.gt.0.0d0 .and. (jp.ne.i+4 .or. jp1.ne.i+3
4439 5097c6509,6511
4440 <                   ecorr6=ecorr6+eello6(i,j,i+1,j1,jj,kk)
4441 ---
4442 >                   ecorr6=ecorr6+eello6(i,jp,i+1,jp1,jj,kk)
4443 >                   if (energy_dec) write (iout,'(a6,4i5,0pf7.3)')
4444 >      1                'ecorr6',i,j,i+1,j1,eello6(i,jp,i+1,jp1,jj,kk)
4445 5101,5103c6515,6517
4446 < cd     &          dabs(eello4(i,j,i+1,j1,jj,kk)),
4447 < cd     &          dabs(eello5(i,j,i+1,j1,jj,kk)),
4448 < cd     &          dabs(eello6(i,j,i+1,j1,jj,kk))
4449 ---
4450 > cd     &          dabs(eello4(i,jp,i+1,jp1,jj,kk)),
4451 > cd     &          dabs(eello5(i,jp,i+1,jp1,jj,kk)),
4452 > cd     &          dabs(eello6(i,jp,i+1,jp1,jj,kk))
4453 5105,5106c6519,6520
4454 <      &            .and. (j.eq.i+4 .and. j1.eq.i+3)) then
4455 < cd                  write (iout,*) '******eturn6: i,j,i+1,j1',i,j,i+1,j1
4456 ---
4457 >      &            .and. (jp.eq.i+4 .and. jp1.eq.i+3)) then
4458 > cd                  write (iout,*) '******eturn6: i,j,i+1,j1',i,jip,i+1,jp1
4459 5107a6522,6523
4460 >                   if (energy_dec) write (iout,'(a6,4i5,0pf7.3)')
4461 >      1                 'eturn6',i,j,i+1,j1,eello_turn6(i,jj,kk)
4462 5112,5115d6527
4463 <             else if (j1.eq.j) then
4464 < C Contacts I-J and I-(J+1) occur simultaneously.
4465 < C The system loses extra energy.
4466 < c             ecorr=ecorr+ehbcorr(i,j,i+1,j,jj,kk,0.60D0,-0.40D0)
4467 5118,5127d6529
4468 <           do kk=1,num_conti
4469 <             j1=jcont_hb(kk,i)
4470 < c           write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1,
4471 < c    &         ' jj=',jj,' kk=',kk
4472 <             if (j1.eq.j+1) then
4473 < C Contacts I-J and (I+1)-J occur simultaneously.
4474 < C The system loses extra energy.
4475 < c             ecorr=ecorr+ehbcorr(i,j,i,j+1,jj,kk,0.60D0,-0.40D0)
4476 <             endif ! j1==j+1
4477 <           enddo ! kk
4478 5129a6532,6594
4479 >       do i=1,nres
4480 >         num_cont_hb(i)=num_cont_hb_old(i)
4481 >       enddo
4482 > c                write (iout,*) "gradcorr5 in eello5"
4483 > c                do iii=1,nres
4484 > c                  write (iout,'(i5,3f10.5)')
4485 > c     &             iii,(gradcorr5(jjj,iii),jjj=1,3)
4486 > c                enddo
4487 >       return
4488 >       end
4489 > c------------------------------------------------------------------------------
4490 >       subroutine add_hb_contact_eello(ii,jj,itask)
4491 >       implicit real*8 (a-h,o-z)
4492 >       include "DIMENSIONS"
4493 >       include "COMMON.IOUNITS"
4494 >       integer max_cont
4495 >       integer max_dim
4496 >       parameter (max_cont=maxconts)
4497 >       parameter (max_dim=70)
4498 >       include "COMMON.CONTACTS"
4499 >       double precision zapas(max_dim,maxconts,max_fg_procs),
4500 >      &  zapas_recv(max_dim,maxconts,max_fg_procs)
4501 >       common /przechowalnia/ zapas
4502 >       integer i,j,ii,jj,iproc,itask(4),nn
4503 > c      write (iout,*) "itask",itask
4504 >       do i=1,2
4505 >         iproc=itask(i)
4506 >         if (iproc.gt.0) then
4507 >           do j=1,num_cont_hb(ii)
4508 >             jjc=jcont_hb(j,ii)
4509 > c            write (iout,*) "send turns i",ii," j",jj," jjc",jjc
4510 >             if (jjc.eq.jj) then
4511 >               ncont_sent(iproc)=ncont_sent(iproc)+1
4512 >               nn=ncont_sent(iproc)
4513 >               zapas(1,nn,iproc)=ii
4514 >               zapas(2,nn,iproc)=jjc
4515 >               zapas(3,nn,iproc)=d_cont(j,ii)
4516 >               ind=3
4517 >               do kk=1,3
4518 >                 ind=ind+1
4519 >                 zapas(ind,nn,iproc)=grij_hb_cont(kk,j,ii)
4520 >               enddo
4521 >               do kk=1,2
4522 >                 do ll=1,2
4523 >                   ind=ind+1
4524 >                   zapas(ind,nn,iproc)=a_chuj(ll,kk,j,ii)
4525 >                 enddo
4526 >               enddo
4527 >               do jj=1,5
4528 >                 do kk=1,3
4529 >                   do ll=1,2
4530 >                     do mm=1,2
4531 >                       ind=ind+1
4532 >                       zapas(ind,nn,iproc)=a_chuj_der(mm,ll,kk,jj,j,ii)
4533 >                     enddo
4534 >                   enddo
4535 >                 enddo
4536 >               enddo
4537 >               exit
4538 >             endif
4539 >           enddo
4540 >         endif
4541 >       enddo
4542 5157,5161c6622,6626
4543 < c     write (iout,*)'Contacts have occurred for peptide groups',i,j,
4544 < c    &   ' and',k,l
4545 < c     write (iout,*)'Contacts have occurred for peptide groups',
4546 < c    &  i,j,' fcont:',eij,' eij',' eesij',ees0pij,ees0mij,' and ',k,l
4547 < c    & ,' fcont ',ekl,' eeskl',ees0pkl,ees0mkl,' ees=',ees
4548 ---
4549 > c      write (iout,'(2(a,2i3,a,f10.5,a,2f10.5),a,f10.5,a,$)')
4550 > c     & 'Contacts ',i,j,
4551 > c     & ' eij',eij,' eesij',ees0pij,ees0mij,' and ',k,l
4552 > c     & ,' fcont ',ekl,' eeskl',ees0pkl,ees0mkl,' energy=',ekont*ees,
4553 > c     & 'gradcorr_long'
4554 5163,5164c6628
4555 <       ecorr=ecorr+ekont*ees
4556 <       if (calc_grad) then
4557 ---
4558 > c      ecorr=ecorr+ekont*ees
4559 5165a6630,6633
4560 >       coeffpees0pij=coeffp*ees0pij
4561 >       coeffmees0mij=coeffm*ees0mij
4562 >       coeffpees0pkl=coeffp*ees0pkl
4563 >       coeffmees0mkl=coeffm*ees0mkl
4564 5167,5198c6635,6678
4565 <         ghalf=0.5D0*ees*ekl*gacont_hbr(ll,jj,i)
4566 <         gradcorr(ll,i)=gradcorr(ll,i)+ghalf
4567 <      &  -ekont*(coeffp*ees0pkl*gacontp_hb1(ll,jj,i)+
4568 <      &  coeffm*ees0mkl*gacontm_hb1(ll,jj,i))
4569 <         gradcorr(ll,j)=gradcorr(ll,j)+ghalf
4570 <      &  -ekont*(coeffp*ees0pkl*gacontp_hb2(ll,jj,i)+
4571 <      &  coeffm*ees0mkl*gacontm_hb2(ll,jj,i))
4572 <         ghalf=0.5D0*ees*eij*gacont_hbr(ll,kk,k)
4573 <         gradcorr(ll,k)=gradcorr(ll,k)+ghalf
4574 <      &  -ekont*(coeffp*ees0pij*gacontp_hb1(ll,kk,k)+
4575 <      &  coeffm*ees0mij*gacontm_hb1(ll,kk,k))
4576 <         gradcorr(ll,l)=gradcorr(ll,l)+ghalf
4577 <      &  -ekont*(coeffp*ees0pij*gacontp_hb2(ll,kk,k)+
4578 <      &  coeffm*ees0mij*gacontm_hb2(ll,kk,k))
4579 <       enddo
4580 <       do m=i+1,j-1
4581 <         do ll=1,3
4582 <           gradcorr(ll,m)=gradcorr(ll,m)+
4583 <      &     ees*ekl*gacont_hbr(ll,jj,i)-
4584 <      &     ekont*(coeffp*ees0pkl*gacontp_hb3(ll,jj,i)+
4585 <      &     coeffm*ees0mkl*gacontm_hb3(ll,jj,i))
4586 <         enddo
4587 <       enddo
4588 <       do m=k+1,l-1
4589 <         do ll=1,3
4590 <           gradcorr(ll,m)=gradcorr(ll,m)+
4591 <      &     ees*eij*gacont_hbr(ll,kk,k)-
4592 <      &     ekont*(coeffp*ees0pij*gacontp_hb3(ll,kk,k)+
4593 <      &     coeffm*ees0mij*gacontm_hb3(ll,kk,k))
4594 <         enddo
4595 <       enddo
4596 <       endif
4597 ---
4598 > cgrad        ghalfi=ees*ekl*gacont_hbr(ll,jj,i)
4599 >         gradcorr(ll,i)=gradcorr(ll,i)!+0.5d0*ghalfi
4600 >      &  -ekont*(coeffpees0pkl*gacontp_hb1(ll,jj,i)+
4601 >      &  coeffmees0mkl*gacontm_hb1(ll,jj,i))
4602 >         gradcorr(ll,j)=gradcorr(ll,j)!+0.5d0*ghalfi
4603 >      &  -ekont*(coeffpees0pkl*gacontp_hb2(ll,jj,i)+
4604 >      &  coeffmees0mkl*gacontm_hb2(ll,jj,i))
4605 > cgrad        ghalfk=ees*eij*gacont_hbr(ll,kk,k)
4606 >         gradcorr(ll,k)=gradcorr(ll,k)!+0.5d0*ghalfk
4607 >      &  -ekont*(coeffpees0pij*gacontp_hb1(ll,kk,k)+
4608 >      &  coeffmees0mij*gacontm_hb1(ll,kk,k))
4609 >         gradcorr(ll,l)=gradcorr(ll,l)!+0.5d0*ghalfk
4610 >      &  -ekont*(coeffpees0pij*gacontp_hb2(ll,kk,k)+
4611 >      &  coeffmees0mij*gacontm_hb2(ll,kk,k))
4612 >         gradlongij=ees*ekl*gacont_hbr(ll,jj,i)-
4613 >      &     ekont*(coeffpees0pkl*gacontp_hb3(ll,jj,i)+
4614 >      &     coeffmees0mkl*gacontm_hb3(ll,jj,i))
4615 >         gradcorr_long(ll,j)=gradcorr_long(ll,j)+gradlongij
4616 >         gradcorr_long(ll,i)=gradcorr_long(ll,i)-gradlongij
4617 >         gradlongkl=ees*eij*gacont_hbr(ll,kk,k)-
4618 >      &     ekont*(coeffpees0pij*gacontp_hb3(ll,kk,k)+
4619 >      &     coeffmees0mij*gacontm_hb3(ll,kk,k))
4620 >         gradcorr_long(ll,l)=gradcorr_long(ll,l)+gradlongkl
4621 >         gradcorr_long(ll,k)=gradcorr_long(ll,k)-gradlongkl
4622 > c        write (iout,'(2f10.5,2x,$)') gradlongij,gradlongkl
4623 >       enddo
4624 > c      write (iout,*)
4625 > cgrad      do m=i+1,j-1
4626 > cgrad        do ll=1,3
4627 > cgrad          gradcorr(ll,m)=gradcorr(ll,m)+
4628 > cgrad     &     ees*ekl*gacont_hbr(ll,jj,i)-
4629 > cgrad     &     ekont*(coeffp*ees0pkl*gacontp_hb3(ll,jj,i)+
4630 > cgrad     &     coeffm*ees0mkl*gacontm_hb3(ll,jj,i))
4631 > cgrad        enddo
4632 > cgrad      enddo
4633 > cgrad      do m=k+1,l-1
4634 > cgrad        do ll=1,3
4635 > cgrad          gradcorr(ll,m)=gradcorr(ll,m)+
4636 > cgrad     &     ees*eij*gacont_hbr(ll,kk,k)-
4637 > cgrad     &     ekont*(coeffp*ees0pij*gacontp_hb3(ll,kk,k)+
4638 > cgrad     &     coeffm*ees0mij*gacontm_hb3(ll,kk,k))
4639 > cgrad        enddo
4640 > cgrad      enddo
4641 > c      write (iout,*) "ehbcorr",ekont*ees
4642 5201a6682
4643 > #ifdef MOMENT
4644 5206d6686
4645 <       include 'DIMENSIONS.ZSCOPT'
4646 5240d6719
4647 <       if (.not.calc_grad) return
4648 5264a6744
4649 > #endif
4650 5273d6752
4651 <       include 'DIMENSIONS.ZSCOPT'
4652 5289a6769,6770
4653 > cd      write (iout,*) "a_chujij",((a_chuj(iii,jjj,jj,i),iii=1,2),jjj=1,2)
4654 > cd      write (iout,*) "a_chujkl",((a_chuj(iii,jjj,kk,k),iii=1,2),jjj=1,2)
4655 5650d7130
4656 <       include 'DIMENSIONS.ZSCOPT'
4657 5671d7150
4658 <       if (calc_grad) then
4659 5713,5718c7192,7197
4660 < cold        ghalf=0.5d0*eel4*ekl*gacont_hbr(ll,jj,i)
4661 <         ggg1(ll)=eel4*g_contij(ll,1)
4662 <         ggg2(ll)=eel4*g_contij(ll,2)
4663 <         ghalf=0.5d0*ggg1(ll)
4664 < cd        ghalf=0.0d0
4665 <         gradcorr(ll,i)=gradcorr(ll,i)+ghalf+ekont*derx(ll,2,1)
4666 ---
4667 > cgrad        ggg1(ll)=eel4*g_contij(ll,1)
4668 > cgrad        ggg2(ll)=eel4*g_contij(ll,2)
4669 >         glongij=eel4*g_contij(ll,1)+ekont*derx(ll,1,1)
4670 >         glongkl=eel4*g_contij(ll,2)+ekont*derx(ll,1,2)
4671 > cgrad        ghalf=0.5d0*ggg1(ll)
4672 >         gradcorr(ll,i)=gradcorr(ll,i)+ekont*derx(ll,2,1)
4673 5720c7199
4674 <         gradcorr(ll,j)=gradcorr(ll,j)+ghalf+ekont*derx(ll,4,1)
4675 ---
4676 >         gradcorr(ll,j)=gradcorr(ll,j)+ekont*derx(ll,4,1)
4677 5722,5725c7201,7204
4678 < cold        ghalf=0.5d0*eel4*eij*gacont_hbr(ll,kk,k)
4679 <         ghalf=0.5d0*ggg2(ll)
4680 < cd        ghalf=0.0d0
4681 <         gradcorr(ll,k)=gradcorr(ll,k)+ghalf+ekont*derx(ll,2,2)
4682 ---
4683 >         gradcorr_long(ll,j)=gradcorr_long(ll,j)+glongij
4684 >         gradcorr_long(ll,i)=gradcorr_long(ll,i)-glongij
4685 > cgrad        ghalf=0.5d0*ggg2(ll)
4686 >         gradcorr(ll,k)=gradcorr(ll,k)+ekont*derx(ll,2,2)
4687 5727c7206
4688 <         gradcorr(ll,l)=gradcorr(ll,l)+ghalf+ekont*derx(ll,4,2)
4689 ---
4690 >         gradcorr(ll,l)=gradcorr(ll,l)+ekont*derx(ll,4,2)
4691 5728a7208,7209
4692 >         gradcorr_long(ll,l)=gradcorr_long(ll,l)+glongkl
4693 >         gradcorr_long(ll,k)=gradcorr_long(ll,k)-glongkl
4694 5730,5753c7211,7230
4695 < cd      goto 1112
4696 <       do m=i+1,j-1
4697 <         do ll=1,3
4698 < cold          gradcorr(ll,m)=gradcorr(ll,m)+eel4*ekl*gacont_hbr(ll,jj,i)
4699 <           gradcorr(ll,m)=gradcorr(ll,m)+ggg1(ll)
4700 <         enddo
4701 <       enddo
4702 <       do m=k+1,l-1
4703 <         do ll=1,3
4704 < cold          gradcorr(ll,m)=gradcorr(ll,m)+eel4*eij*gacont_hbr(ll,kk,k)
4705 <           gradcorr(ll,m)=gradcorr(ll,m)+ggg2(ll)
4706 <         enddo
4707 <       enddo
4708 < 1112  continue
4709 <       do m=i+2,j2
4710 <         do ll=1,3
4711 <           gradcorr(ll,m)=gradcorr(ll,m)+ekont*derx(ll,1,1)
4712 <         enddo
4713 <       enddo
4714 <       do m=k+2,l2
4715 <         do ll=1,3
4716 <           gradcorr(ll,m)=gradcorr(ll,m)+ekont*derx(ll,1,2)
4717 <         enddo
4718 <       enddo
4719 ---
4720 > cgrad      do m=i+1,j-1
4721 > cgrad        do ll=1,3
4722 > cgrad          gradcorr(ll,m)=gradcorr(ll,m)+ggg1(ll)
4723 > cgrad        enddo
4724 > cgrad      enddo
4725 > cgrad      do m=k+1,l-1
4726 > cgrad        do ll=1,3
4727 > cgrad          gradcorr(ll,m)=gradcorr(ll,m)+ggg2(ll)
4728 > cgrad        enddo
4729 > cgrad      enddo
4730 > cgrad      do m=i+2,j2
4731 > cgrad        do ll=1,3
4732 > cgrad          gradcorr(ll,m)=gradcorr(ll,m)+ekont*derx(ll,1,1)
4733 > cgrad        enddo
4734 > cgrad      enddo
4735 > cgrad      do m=k+2,l2
4736 > cgrad        do ll=1,3
4737 > cgrad          gradcorr(ll,m)=gradcorr(ll,m)+ekont*derx(ll,1,2)
4738 > cgrad        enddo
4739 > cgrad      enddo
4740 5757d7233
4741 <       endif
4742 5767d7242
4743 <       include 'DIMENSIONS.ZSCOPT'
4744 5847d7321
4745 <       if (calc_grad) then
4746 5886d7359
4747 <       endif
4748 5895d7367
4749 <       if (calc_grad) then
4750 5926d7397
4751 <       endif
4752 5938d7408
4753 <         if (calc_grad) then
4754 5971d7440
4755 <         endif
4756 5980d7448
4757 <         if (calc_grad) then
4758 6004d7471
4759 <         endif
4760 6015d7481
4761 <         if (calc_grad) then
4762 6048d7513
4763 <         endif
4764 6057d7521
4765 <         if (calc_grad) then
4766 6082d7545
4767 <       endif
4768 6094d7556
4769 <       if (calc_grad) then
4770 6112a7575,7578
4771 > C 2/11/08 AL Gradients over DC's connecting interacting sites will be
4772 > C        summed up outside the subrouine as for the other subroutines
4773 > C        handling long-range interactions. The old code is commented out
4774 > C        with "cgrad" to keep track of changes.
4775 6114,6115c7580,7591
4776 <         ggg1(ll)=eel5*g_contij(ll,1)
4777 <         ggg2(ll)=eel5*g_contij(ll,2)
4778 ---
4779 > cgrad        ggg1(ll)=eel5*g_contij(ll,1)
4780 > cgrad        ggg2(ll)=eel5*g_contij(ll,2)
4781 >         gradcorr5ij=eel5*g_contij(ll,1)+ekont*derx(ll,1,1)
4782 >         gradcorr5kl=eel5*g_contij(ll,2)+ekont*derx(ll,1,2)
4783 > c        write (iout,'(a,3i3,a,5f8.3,2i3,a,5f8.3,a,f8.3)')
4784 > c     &   "ecorr5",ll,i,j," derx",derx(ll,2,1),derx(ll,3,1),derx(ll,4,1),
4785 > c     &   derx(ll,5,1),k,l," derx",derx(ll,2,2),derx(ll,3,2),
4786 > c     &   derx(ll,4,2),derx(ll,5,2)," ekont",ekont
4787 > c        write (iout,'(a,3i3,a,3f8.3,2i3,a,3f8.3)')
4788 > c     &   "ecorr5",ll,i,j," gradcorr5",g_contij(ll,1),derx(ll,1,1),
4789 > c     &   gradcorr5ij,
4790 > c     &   k,l," gradcorr5",g_contij(ll,2),derx(ll,1,2),gradcorr5kl
4791 6117c7593
4792 <         ghalf=0.5d0*ggg1(ll)
4793 ---
4794 > cgrad        ghalf=0.5d0*ggg1(ll)
4795 6119c7595
4796 <         gradcorr5(ll,i)=gradcorr5(ll,i)+ghalf+ekont*derx(ll,2,1)
4797 ---
4798 >         gradcorr5(ll,i)=gradcorr5(ll,i)+ekont*derx(ll,2,1)
4799 6121c7597
4800 <         gradcorr5(ll,j)=gradcorr5(ll,j)+ghalf+ekont*derx(ll,4,1)
4801 ---
4802 >         gradcorr5(ll,j)=gradcorr5(ll,j)+ekont*derx(ll,4,1)
4803 6122a7599,7600
4804 >         gradcorr5_long(ll,j)=gradcorr5_long(ll,j)+gradcorr5ij
4805 >         gradcorr5_long(ll,i)=gradcorr5_long(ll,i)-gradcorr5ij
4806 6124c7602
4807 <         ghalf=0.5d0*ggg2(ll)
4808 ---
4809 > cgrad        ghalf=0.5d0*ggg2(ll)
4810 6129a7608,7609
4811 >         gradcorr5_long(ll,l)=gradcorr5_long(ll,l)+gradcorr5kl
4812 >         gradcorr5_long(ll,k)=gradcorr5_long(ll,k)-gradcorr5kl
4813 6132,6133c7612,7613
4814 <       do m=i+1,j-1
4815 <         do ll=1,3
4816 ---
4817 > cgrad      do m=i+1,j-1
4818 > cgrad        do ll=1,3
4819 6135,6139c7615,7619
4820 <           gradcorr5(ll,m)=gradcorr5(ll,m)+ggg1(ll)
4821 <         enddo
4822 <       enddo
4823 <       do m=k+1,l-1
4824 <         do ll=1,3
4825 ---
4826 > cgrad          gradcorr5(ll,m)=gradcorr5(ll,m)+ggg1(ll)
4827 > cgrad        enddo
4828 > cgrad      enddo
4829 > cgrad      do m=k+1,l-1
4830 > cgrad        do ll=1,3
4831 6141,6143c7621,7623
4832 <           gradcorr5(ll,m)=gradcorr5(ll,m)+ggg2(ll)
4833 <         enddo
4834 <       enddo
4835 ---
4836 > cgrad          gradcorr5(ll,m)=gradcorr5(ll,m)+ggg2(ll)
4837 > cgrad        enddo
4838 > cgrad      enddo
4839 6145,6154c7625,7634
4840 <       do m=i+2,j2
4841 <         do ll=1,3
4842 <           gradcorr5(ll,m)=gradcorr5(ll,m)+ekont*derx(ll,1,1)
4843 <         enddo
4844 <       enddo
4845 <       do m=k+2,l2
4846 <         do ll=1,3
4847 <           gradcorr5(ll,m)=gradcorr5(ll,m)+ekont*derx(ll,1,2)
4848 <         enddo
4849 <       enddo
4850 ---
4851 > cgrad      do m=i+2,j2
4852 > cgrad        do ll=1,3
4853 > cgrad          gradcorr5(ll,m)=gradcorr5(ll,m)+ekont*derx(ll,1,1)
4854 > cgrad        enddo
4855 > cgrad      enddo
4856 > cgrad      do m=k+2,l2
4857 > cgrad        do ll=1,3
4858 > cgrad          gradcorr5(ll,m)=gradcorr5(ll,m)+ekont*derx(ll,1,2)
4859 > cgrad        enddo
4860 > cgrad      enddo
4861 6158d7637
4862 <       endif
4863 6168d7646
4864 <       include 'DIMENSIONS.ZSCOPT'
4865 6228,6233c7706,7711
4866 < cd      write(iout,*) 'eello6_1',eello6_1,' eel6_1_num',16*eel6_1_num
4867 < cd      write(iout,*) 'eello6_2',eello6_2,' eel6_2_num',16*eel6_2_num
4868 < cd      write(iout,*) 'eello6_3',eello6_3,' eel6_3_num',16*eel6_3_num
4869 < cd      write(iout,*) 'eello6_4',eello6_4,' eel6_4_num',16*eel6_4_num
4870 < cd      write(iout,*) 'eello6_5',eello6_5,' eel6_5_num',16*eel6_5_num
4871 < cd      write(iout,*) 'eello6_6',eello6_6,' eel6_6_num',16*eel6_6_num
4872 ---
4873 > cd      write(iout,*) 'eello6_1',eello6_1!,' eel6_1_num',16*eel6_1_num
4874 > cd      write(iout,*) 'eello6_2',eello6_2!,' eel6_2_num',16*eel6_2_num
4875 > cd      write(iout,*) 'eello6_3',eello6_3!,' eel6_3_num',16*eel6_3_num
4876 > cd      write(iout,*) 'eello6_4',eello6_4!,' eel6_4_num',16*eel6_4_num
4877 > cd      write(iout,*) 'eello6_5',eello6_5!,' eel6_5_num',16*eel6_5_num
4878 > cd      write(iout,*) 'eello6_6',eello6_6!,' eel6_6_num',16*eel6_6_num
4879 6235d7712
4880 <       if (calc_grad) then
4881 6251,6252c7728,7729
4882 <         ggg1(ll)=eel6*g_contij(ll,1)
4883 <         ggg2(ll)=eel6*g_contij(ll,2)
4884 ---
4885 > cgrad        ggg1(ll)=eel6*g_contij(ll,1)
4886 > cgrad        ggg2(ll)=eel6*g_contij(ll,2)
4887 6254c7731
4888 <         ghalf=0.5d0*ggg1(ll)
4889 ---
4890 > cgrad        ghalf=0.5d0*ggg1(ll)
4891 6256c7733,7735
4892 <         gradcorr6(ll,i)=gradcorr6(ll,i)+ghalf+ekont*derx(ll,2,1)
4893 ---
4894 >         gradcorr6ij=eel6*g_contij(ll,1)+ekont*derx(ll,1,1)
4895 >         gradcorr6kl=eel6*g_contij(ll,2)+ekont*derx(ll,1,2)
4896 >         gradcorr6(ll,i)=gradcorr6(ll,i)+ekont*derx(ll,2,1)
4897 6258c7737
4898 <         gradcorr6(ll,j)=gradcorr6(ll,j)+ghalf+ekont*derx(ll,4,1)
4899 ---
4900 >         gradcorr6(ll,j)=gradcorr6(ll,j)+ekont*derx(ll,4,1)
4901 6260c7739,7741
4902 <         ghalf=0.5d0*ggg2(ll)
4903 ---
4904 >         gradcorr6_long(ll,j)=gradcorr6_long(ll,j)+gradcorr6ij
4905 >         gradcorr6_long(ll,i)=gradcorr6_long(ll,i)-gradcorr6ij
4906 > cgrad        ghalf=0.5d0*ggg2(ll)
4907 6263c7744
4908 <         gradcorr6(ll,k)=gradcorr6(ll,k)+ghalf+ekont*derx(ll,2,2)
4909 ---
4910 >         gradcorr6(ll,k)=gradcorr6(ll,k)+ekont*derx(ll,2,2)
4911 6265c7746
4912 <         gradcorr6(ll,l)=gradcorr6(ll,l)+ghalf+ekont*derx(ll,4,2)
4913 ---
4914 >         gradcorr6(ll,l)=gradcorr6(ll,l)+ekont*derx(ll,4,2)
4915 6266a7748,7749
4916 >         gradcorr6_long(ll,l)=gradcorr6_long(ll,l)+gradcorr6kl
4917 >         gradcorr6_long(ll,k)=gradcorr6_long(ll,k)-gradcorr6kl
4918 6269,6270c7752,7753
4919 <       do m=i+1,j-1
4920 <         do ll=1,3
4921 ---
4922 > cgrad      do m=i+1,j-1
4923 > cgrad        do ll=1,3
4924 6272,6276c7755,7759
4925 <           gradcorr6(ll,m)=gradcorr6(ll,m)+ggg1(ll)
4926 <         enddo
4927 <       enddo
4928 <       do m=k+1,l-1
4929 <         do ll=1,3
4930 ---
4931 > cgrad          gradcorr6(ll,m)=gradcorr6(ll,m)+ggg1(ll)
4932 > cgrad        enddo
4933 > cgrad      enddo
4934 > cgrad      do m=k+1,l-1
4935 > cgrad        do ll=1,3
4936 6278,6291c7761,7774
4937 <           gradcorr6(ll,m)=gradcorr6(ll,m)+ggg2(ll)
4938 <         enddo
4939 <       enddo
4940 < 1112  continue
4941 <       do m=i+2,j2
4942 <         do ll=1,3
4943 <           gradcorr6(ll,m)=gradcorr6(ll,m)+ekont*derx(ll,1,1)
4944 <         enddo
4945 <       enddo
4946 <       do m=k+2,l2
4947 <         do ll=1,3
4948 <           gradcorr6(ll,m)=gradcorr6(ll,m)+ekont*derx(ll,1,2)
4949 <         enddo
4950 <       enddo
4951 ---
4952 > cgrad          gradcorr6(ll,m)=gradcorr6(ll,m)+ggg2(ll)
4953 > cgrad        enddo
4954 > cgrad      enddo
4955 > cgrad1112  continue
4956 > cgrad      do m=i+2,j2
4957 > cgrad        do ll=1,3
4958 > cgrad          gradcorr6(ll,m)=gradcorr6(ll,m)+ekont*derx(ll,1,1)
4959 > cgrad        enddo
4960 > cgrad      enddo
4961 > cgrad      do m=k+2,l2
4962 > cgrad        do ll=1,3
4963 > cgrad          gradcorr6(ll,m)=gradcorr6(ll,m)+ekont*derx(ll,1,2)
4964 > cgrad        enddo
4965 > cgrad      enddo
4966 6295d7777
4967 <       endif
4968 6305d7786
4969 <       include 'DIMENSIONS.ZSCOPT'
4970 6346d7826
4971 <       if (.not. calc_grad) return
4972 6411d7890
4973 <       include 'DIMENSIONS.ZSCOPT'
4974 6461d7939
4975 <       if (.not. calc_grad) return
4976 6598d8075
4977 <       include 'DIMENSIONS.ZSCOPT'
4978 6651c8128,8129
4979 < cd      write (2,*) 'eello6_graph3:','s1',s1,' s2',s2,' s3',s3,' s4',s4
4980 ---
4981 > cd      write (2,*) 'eello6_graph3:','s1',s1,' s2',s2,' s3',s3,' s4',s4,
4982 > cd     & "sum",-(s2+s3+s4)
4983 6658d8135
4984 <       if (.not. calc_grad) return
4985 6714d8190
4986 <       include 'DIMENSIONS.ZSCOPT'
4987 6794d8269
4988 <       if (.not. calc_grad) return
4989 6960d8434
4990 <       include 'DIMENSIONS.ZSCOPT'
4991 6975a8450,8453
4992 >       s1=0.0d0
4993 >       s8=0.0d0
4994 >       s13=0.0d0
4995 > c
4996 7013,7014d8490
4997 < #else
4998 <       s1 = 0.0d0
4999 7024,7025d8499
5000 < #else
5001 <       s8=0.0d0
5002 7037,7038d8510
5003 < #else
5004 <       s13=0.0d0
5005 7047d8518
5006 <       if (calc_grad) then
5007 7048a8520,8521
5008 >       s1d =0.0d0
5009 >       s8d =0.0d0
5010 7056,7057d8528
5011 < #else
5012 <       s8d=0.0d0
5013 7074,7075d8544
5014 < #else
5015 <       s1d=0.0d0
5016 7089,7090d8557
5017 < #else
5018 <       s13d=0.0d0
5019 7112,7113d8578
5020 < #else
5021 <       s13d = 0.0d0
5022 7130,7131d8594
5023 < #else
5024 <       s1d = 0.0d0
5025 7140,7141d8602
5026 < #else
5027 <       s8d = 0.0d0
5028 7149,7150d8609
5029 < #else
5030 <       s13d = 0.0d0
5031 7172,7173d8630
5032 < #else
5033 <             s1d = 0.0d0
5034 7184,7185d8640
5035 < #else
5036 <             s8d = 0.0d0
5037 7248,7250c8703,8705
5038 <         ggg1(ll)=eel_turn6*g_contij(ll,1)
5039 <         ggg2(ll)=eel_turn6*g_contij(ll,2)
5040 <         ghalf=0.5d0*ggg1(ll)
5041 ---
5042 > cgrad        ggg1(ll)=eel_turn6*g_contij(ll,1)
5043 > cgrad        ggg2(ll)=eel_turn6*g_contij(ll,2)
5044 > cgrad        ghalf=0.5d0*ggg1(ll)
5045 7252c8707,8709
5046 <         gcorr6_turn(ll,i)=gcorr6_turn(ll,i)+ghalf
5047 ---
5048 >         gturn6ij=eel_turn6*g_contij(ll,1)+ekont*derx_turn(ll,1,1)
5049 >         gturn6kl=eel_turn6*g_contij(ll,2)+ekont*derx_turn(ll,1,2)
5050 >         gcorr6_turn(ll,i)=gcorr6_turn(ll,i)!+ghalf
5051 7255c8712
5052 <         gcorr6_turn(ll,j)=gcorr6_turn(ll,j)+ghalf
5053 ---
5054 >         gcorr6_turn(ll,j)=gcorr6_turn(ll,j)!+ghalf
5055 7258c8715,8717
5056 <         ghalf=0.5d0*ggg2(ll)
5057 ---
5058 >         gcorr6_turn_long(ll,j)=gcorr6_turn_long(ll,j)+gturn6ij
5059 >         gcorr6_turn_long(ll,i)=gcorr6_turn_long(ll,i)-gturn6ij
5060 > cgrad        ghalf=0.5d0*ggg2(ll)
5061 7260c8719
5062 <         gcorr6_turn(ll,k)=gcorr6_turn(ll,k)+ghalf
5063 ---
5064 >         gcorr6_turn(ll,k)=gcorr6_turn(ll,k)!+ghalf
5065 7263c8722
5066 <         gcorr6_turn(ll,l)=gcorr6_turn(ll,l)+ghalf
5067 ---
5068 >         gcorr6_turn(ll,l)=gcorr6_turn(ll,l)!+ghalf
5069 7265a8725,8726
5070 >         gcorr6_turn_long(ll,l)=gcorr6_turn_long(ll,l)+gturn6kl
5071 >         gcorr6_turn_long(ll,k)=gcorr6_turn_long(ll,k)-gturn6kl
5072 7268,7288c8729,8749
5073 <       do m=i+1,j-1
5074 <         do ll=1,3
5075 <           gcorr6_turn(ll,m)=gcorr6_turn(ll,m)+ggg1(ll)
5076 <         enddo
5077 <       enddo
5078 <       do m=k+1,l-1
5079 <         do ll=1,3
5080 <           gcorr6_turn(ll,m)=gcorr6_turn(ll,m)+ggg2(ll)
5081 <         enddo
5082 <       enddo
5083 < 1112  continue
5084 <       do m=i+2,j2
5085 <         do ll=1,3
5086 <           gcorr6_turn(ll,m)=gcorr6_turn(ll,m)+ekont*derx_turn(ll,1,1)
5087 <         enddo
5088 <       enddo
5089 <       do m=k+2,l2
5090 <         do ll=1,3
5091 <           gcorr6_turn(ll,m)=gcorr6_turn(ll,m)+ekont*derx_turn(ll,1,2)
5092 <         enddo
5093 <       enddo
5094 ---
5095 > cgrad      do m=i+1,j-1
5096 > cgrad        do ll=1,3
5097 > cgrad          gcorr6_turn(ll,m)=gcorr6_turn(ll,m)+ggg1(ll)
5098 > cgrad        enddo
5099 > cgrad      enddo
5100 > cgrad      do m=k+1,l-1
5101 > cgrad        do ll=1,3
5102 > cgrad          gcorr6_turn(ll,m)=gcorr6_turn(ll,m)+ggg2(ll)
5103 > cgrad        enddo
5104 > cgrad      enddo
5105 > cgrad1112  continue
5106 > cgrad      do m=i+2,j2
5107 > cgrad        do ll=1,3
5108 > cgrad          gcorr6_turn(ll,m)=gcorr6_turn(ll,m)+ekont*derx_turn(ll,1,1)
5109 > cgrad        enddo
5110 > cgrad      enddo
5111 > cgrad      do m=k+2,l2
5112 > cgrad        do ll=1,3
5113 > cgrad          gcorr6_turn(ll,m)=gcorr6_turn(ll,m)+ekont*derx_turn(ll,1,2)
5114 > cgrad        enddo
5115 > cgrad      enddo
5116 7292d8752
5117 <       endif
5118 7297a8758,8777
5119 >
5120 > C-----------------------------------------------------------------------------
5121 >       double precision function scalar(u,v)
5122 > !DIR$ INLINEALWAYS scalar
5123 > #ifndef OSF
5124 > cDEC$ ATTRIBUTES FORCEINLINE::scalar
5125 > #endif
5126 >       implicit none
5127 >       double precision u(3),v(3)
5128 > cd      double precision sc
5129 > cd      integer i
5130 > cd      sc=0.0d0
5131 > cd      do i=1,3
5132 > cd        sc=sc+u(i)*v(i)
5133 > cd      enddo
5134 > cd      scalar=sc
5135 >
5136 >       scalar=u(1)*v(1)+u(2)*v(2)+u(3)*v(3)
5137 >       return
5138 >       end
5139 7299a8780,8783
5140 > !DIR$ INLINEALWAYS MATVEC2
5141 > #ifndef OSF
5142 > cDEC$ ATTRIBUTES FORCEINLINE::MATVEC2
5143 > #endif
5144 7317a8802,8804
5145 > #ifndef OSF
5146 > cDEC$ ATTRIBUTES FORCEINLINE::MATMAT2
5147 > #endif
5148 7343a8831
5149 > !DIR$ INLINEALWAYS scalar2
5150 7354a8843,8846
5151 > !DIR$ INLINEALWAYS transpose2
5152 > #ifndef OSF
5153 > cDEC$ ATTRIBUTES FORCEINLINE::transpose2
5154 > #endif
5155 7376a8869,8872
5156 > !DIR$ INLINEALWAYS prodmat3
5157 > #ifndef OSF
5158 > cDEC$ ATTRIBUTES FORCEINLINE::prodmat3
5159 > #endif
5160 7419,7431d8914
5161 < C-----------------------------------------------------------------------------
5162 <       double precision function scalar(u,v)
5163 <       implicit none
5164 <       double precision u(3),v(3)
5165 <       double precision sc
5166 <       integer i
5167 <       sc=0.0d0
5168 <       do i=1,3
5169 <         sc=sc+u(i)*v(i)
5170 <       enddo
5171 <       scalar=sc
5172 <       return
5173 <       end