source/unres/src_MD-M/sc_move.F

   1       subroutine sc_move(n_start,n_end,n_maxtry,e_drop,
   2      +     n_fun,etot)
   3 c     Perform a quick search over side-chain arrangments (over
   4 c     residues n_start to n_end) for a given (frozen) CA trace
   5 c     Only side-chains are minimized (at most n_maxtry times each),
   6 c     not CA positions
   7 c     Stops if energy drops by e_drop, otherwise tries all residues
   8 c     in the given range
   9 c     If there is an energy drop, full minimization may be useful
  10 c     n_start, n_end CAN be modified by this routine, but only if
  11 c     out of bounds (n_start <= 1, n_end >= nres, n_start < n_end)
  12 c     NOTE: this move should never increase the energy
  13 crc      implicit none
  14
  15 c     Includes
  16       implicit real*8 (a-h,o-z)
  17       include 'DIMENSIONS'
  18       include 'mpif.h'
  19       include 'COMMON.GEO'
  20       include 'COMMON.VAR'
  21       include 'COMMON.HEADER'
  22       include 'COMMON.IOUNITS'
  23       include 'COMMON.CHAIN'
  24       include 'COMMON.FFIELD'
  25
  26 c     External functions
  27       integer iran_num
  28       external iran_num
  29
  30 c     Input arguments
  31       integer n_start,n_end,n_maxtry
  32       double precision e_drop
  33
  34 c     Output arguments
  35       integer n_fun
  36       double precision etot
  37
  38 c     Local variables
  39       double precision energy(0:n_ene)
  40       double precision cur_alph(2:nres-1),cur_omeg(2:nres-1)
  41       double precision orig_e,cur_e
  42       integer n,n_steps,n_first,n_cur,n_tot,i
  43       double precision orig_w(n_ene)
  44       double precision wtime
  45
  46
  47 c     Set non side-chain weights to zero (minimization is faster)
  48 c     NOTE: e(2) does not actually depend on the side-chain, only CA
  49       orig_w(2)=wscp
  50       orig_w(3)=welec
  51       orig_w(4)=wcorr
  52       orig_w(5)=wcorr5
  53       orig_w(6)=wcorr6
  54       orig_w(7)=wel_loc
  55       orig_w(8)=wturn3
  56       orig_w(9)=wturn4
  57       orig_w(10)=wturn6
  58       orig_w(11)=wang
  59       orig_w(13)=wtor
  60       orig_w(14)=wtor_d
  61       orig_w(15)=wvdwpp
  62
  63       wscp=0.D0
  64       welec=0.D0
  65       wcorr=0.D0
  66       wcorr5=0.D0
  67       wcorr6=0.D0
  68       wel_loc=0.D0
  69       wturn3=0.D0
  70       wturn4=0.D0
  71       wturn6=0.D0
  72       wang=0.D0
  73       wtor=0.D0
  74       wtor_d=0.D0
  75       wvdwpp=0.D0
  76
  77 c     Make sure n_start, n_end are within proper range
  78       if (n_start.lt.2) n_start=2
  79       if (n_end.gt.nres-1) n_end=nres-1
  80 crc      if (n_start.lt.n_end) then
  81       if (n_start.gt.n_end) then
  82         n_start=2
  83         n_end=nres-1
  84       endif
  85
  86 c     Save the initial values of energy and coordinates
  87 cd      call chainbuild
  88 cd      call etotal(energy)
  89 cd      write (iout,*) 'start sc ene',energy(0)
  90 cd      call enerprint(energy(0))
  91 crc      etot=energy(0)
  92        n_fun=0
  93 crc      orig_e=etot
  94 crc      cur_e=orig_e
  95 crc      do i=2,nres-1
  96 crc        cur_alph(i)=alph(i)
  97 crc        cur_omeg(i)=omeg(i)
  98 crc      enddo
  99
 100 ct      wtime=MPI_WTIME()
 101 c     Try (one by one) all specified residues, starting from a
 102 c     random position in sequence
 103 c     Stop early if the energy has decreased by at least e_drop
 104       n_tot=n_end-n_start+1
 105       n_first=iran_num(0,n_tot-1)
 106       n_steps=0
 107       n=0
 108 crc      do while (n.lt.n_tot .and. orig_e-etot.lt.e_drop)
 109       do while (n.lt.n_tot)
 110         n_cur=n_start+mod(n_first+n,n_tot)
 111         call single_sc_move(n_cur,n_maxtry,e_drop,
 112      +       n_steps,n_fun,etot)
 113 c     If a lower energy was found, update the current structure...
 114 crc        if (etot.lt.cur_e) then
 115 crc          cur_e=etot
 116 crc          do i=2,nres-1
 117 crc            cur_alph(i)=alph(i)
 118 crc            cur_omeg(i)=omeg(i)
 119 crc          enddo
 120 crc        else
 121 c     ...else revert to the previous one
 122 crc          etot=cur_e
 123 crc          do i=2,nres-1
 124 crc            alph(i)=cur_alph(i)
 125 crc            omeg(i)=cur_omeg(i)
 126 crc          enddo
 127 crc        endif
 128         n=n+1
 129 cd
 130 cd      call chainbuild
 131 cd      call etotal(energy)
 132 cd      print *,'running',n,energy(0)
 133       enddo
 134
 135 cd      call chainbuild
 136 cd      call etotal(energy)
 137 cd      write (iout,*) 'end   sc ene',energy(0)
 138
 139 c     Put the original weights back to calculate the full energy
 140       wscp=orig_w(2)
 141       welec=orig_w(3)
 142       wcorr=orig_w(4)
 143       wcorr5=orig_w(5)
 144       wcorr6=orig_w(6)
 145       wel_loc=orig_w(7)
 146       wturn3=orig_w(8)
 147       wturn4=orig_w(9)
 148       wturn6=orig_w(10)
 149       wang=orig_w(11)
 150       wtor=orig_w(13)
 151       wtor_d=orig_w(14)
 152       wvdwpp=orig_w(15)
 153
 154 crc      n_fun=n_fun+1
 155 ct      write (iout,*) 'sc_local time= ',MPI_WTIME()-wtime
 156       return
 157       end
 158
 159 c-------------------------------------------------------------
 160
 161       subroutine single_sc_move(res_pick,n_maxtry,e_drop,
 162      +     n_steps,n_fun,e_sc)
 163 c     Perturb one side-chain (res_pick) and minimize the
 164 c     neighbouring region, keeping all CA's and non-neighbouring
 165 c     side-chains fixed
 166 c     Try until e_drop energy improvement is achieved, or n_maxtry
 167 c     attempts have been made
 168 c     At the start, e_sc should contain the side-chain-only energy(0)
 169 c     nsteps and nfun for this move are ADDED to n_steps and n_fun
 170 crc      implicit none
 171
 172 c     Includes
 173       implicit real*8 (a-h,o-z)
 174       include 'DIMENSIONS'
 175       include 'COMMON.VAR'
 176       include 'COMMON.INTERACT'
 177       include 'COMMON.CHAIN'
 178       include 'COMMON.MINIM'
 179       include 'COMMON.FFIELD'
 180       include 'COMMON.IOUNITS'
 181
 182 c     External functions
 183       double precision dist
 184       external dist
 185
 186 c     Input arguments
 187       integer res_pick,n_maxtry
 188       double precision e_drop
 189
 190 c     Input/Output arguments
 191       integer n_steps,n_fun
 192       double precision e_sc
 193
 194 c     Local variables
 195       logical fail
 196       integer i,j
 197       integer nres_moved
 198       integer iretcode,loc_nfun,orig_maxfun,n_try
 199       double precision sc_dist,sc_dist_cutoff
 200       double precision energy(0:n_ene),orig_e,cur_e
 201       double precision evdw,escloc
 202       double precision cur_alph(2:nres-1),cur_omeg(2:nres-1)
 203       double precision var(maxvar)
 204
 205       double precision orig_theta(1:nres),orig_phi(1:nres),
 206      +     orig_alph(1:nres),orig_omeg(1:nres)
 207
 208
 209 c     Define what is meant by "neighbouring side-chain"
 210       sc_dist_cutoff=5.0D0
 211
 212 c     Don't do glycine or ends
 213       i=itype(res_pick)
 214       if (i.eq.10 .or. i.eq.21) return
 215
 216 c     Freeze everything (later will relax only selected side-chains)
 217       mask_r=.true.
 218       do i=1,nres
 219         mask_phi(i)=0
 220         mask_theta(i)=0
 221         mask_side(i)=0
 222       enddo
 223
 224 c     Find the neighbours of the side-chain to move
 225 c     and save initial variables
 226 crc      orig_e=e_sc
 227 crc      cur_e=orig_e
 228       nres_moved=0
 229       do i=2,nres-1
 230 c     Don't do glycine (itype(j)==10)
 231         if (itype(i).ne.10) then
 232           sc_dist=dist(nres+i,nres+res_pick)
 233         else
 234           sc_dist=sc_dist_cutoff
 235         endif
 236         if (sc_dist.lt.sc_dist_cutoff) then
 237           nres_moved=nres_moved+1
 238           mask_side(i)=1
 239           cur_alph(i)=alph(i)
 240           cur_omeg(i)=omeg(i)
 241         endif
 242       enddo
 243
 244       call chainbuild
 245       call egb1(evdw)
 246       call esc(escloc)
 247       e_sc=wsc*evdw+wscloc*escloc
 248 cd      call etotal(energy)
 249 cd      print *,'new       ',(energy(k),k=0,n_ene)
 250       orig_e=e_sc
 251       cur_e=orig_e
 252
 253       n_try=0
 254       do while (n_try.lt.n_maxtry .and. orig_e-cur_e.lt.e_drop)
 255 c     Move the selected residue (don't worry if it fails)
 256         call gen_side(itype(res_pick),theta(res_pick+1),
 257      +       alph(res_pick),omeg(res_pick),fail)
 258
 259 c     Minimize the side-chains starting from the new arrangement
 260         call geom_to_var(nvar,var)
 261         orig_maxfun=maxfun
 262         maxfun=7
 263
 264 crc        do i=1,nres
 265 crc          orig_theta(i)=theta(i)
 266 crc          orig_phi(i)=phi(i)
 267 crc          orig_alph(i)=alph(i)
 268 crc          orig_omeg(i)=omeg(i)
 269 crc        enddo
 270
 271         call minimize_sc1(e_sc,var,iretcode,loc_nfun)
 272
 273 cv        write(*,'(2i3,2f12.5,2i3)')
 274 cv     &       res_pick,nres_moved,orig_e,e_sc-cur_e,
 275 cv     &       iretcode,loc_nfun
 276
 277 c$$$        if (iretcode.eq.8) then
 278 c$$$          write(iout,*)'Coordinates just after code 8'
 279 c$$$          call chainbuild
 280 c$$$          call all_varout
 281 c$$$          call flush(iout)
 282 c$$$          do i=1,nres
 283 c$$$            theta(i)=orig_theta(i)
 284 c$$$            phi(i)=orig_phi(i)
 285 c$$$            alph(i)=orig_alph(i)
 286 c$$$            omeg(i)=orig_omeg(i)
 287 c$$$          enddo
 288 c$$$          write(iout,*)'Coordinates just before code 8'
 289 c$$$          call chainbuild
 290 c$$$          call all_varout
 291 c$$$          call flush(iout)
 292 c$$$        endif
 293
 294         n_fun=n_fun+loc_nfun
 295         maxfun=orig_maxfun
 296         call var_to_geom(nvar,var)
 297
 298 c     If a lower energy was found, update the current structure...
 299         if (e_sc.lt.cur_e) then
 300 cv              call chainbuild
 301 cv              call etotal(energy)
 302 cd              call egb1(evdw)
 303 cd              call esc(escloc)
 304 cd              e_sc1=wsc*evdw+wscloc*escloc
 305 cd              print *,'     new',e_sc1,energy(0)
 306 cv              print *,'new       ',energy(0)
 307 cd              call enerprint(energy(0))
 308           cur_e=e_sc
 309           do i=2,nres-1
 310             if (mask_side(i).eq.1) then
 311               cur_alph(i)=alph(i)
 312               cur_omeg(i)=omeg(i)
 313             endif
 314           enddo
 315         else
 316 c     ...else revert to the previous one
 317           e_sc=cur_e
 318           do i=2,nres-1
 319             if (mask_side(i).eq.1) then
 320               alph(i)=cur_alph(i)
 321               omeg(i)=cur_omeg(i)
 322             endif
 323           enddo
 324         endif
 325         n_try=n_try+1
 326
 327       enddo
 328       n_steps=n_steps+n_try
 329
 330 c     Reset the minimization mask_r to false
 331       mask_r=.false.
 332
 333       return
 334       end
 335
 336 c-------------------------------------------------------------
 337
 338       subroutine sc_minimize(etot,iretcode,nfun)
 339 c     Minimizes side-chains only, leaving backbone frozen
 340 crc      implicit none
 341
 342 c     Includes
 343       implicit real*8 (a-h,o-z)
 344       include 'DIMENSIONS'
 345       include 'COMMON.VAR'
 346       include 'COMMON.CHAIN'
 347       include 'COMMON.FFIELD'
 348
 349 c     Output arguments
 350       double precision etot
 351       integer iretcode,nfun
 352
 353 c     Local variables
 354       integer i
 355       double precision orig_w(n_ene),energy(0:n_ene)
 356       double precision var(maxvar)
 357
 358
 359 c     Set non side-chain weights to zero (minimization is faster)
 360 c     NOTE: e(2) does not actually depend on the side-chain, only CA
 361       orig_w(2)=wscp
 362       orig_w(3)=welec
 363       orig_w(4)=wcorr
 364       orig_w(5)=wcorr5
 365       orig_w(6)=wcorr6
 366       orig_w(7)=wel_loc
 367       orig_w(8)=wturn3
 368       orig_w(9)=wturn4
 369       orig_w(10)=wturn6
 370       orig_w(11)=wang
 371       orig_w(13)=wtor
 372       orig_w(14)=wtor_d
 373
 374       wscp=0.D0
 375       welec=0.D0
 376       wcorr=0.D0
 377       wcorr5=0.D0
 378       wcorr6=0.D0
 379       wel_loc=0.D0
 380       wturn3=0.D0
 381       wturn4=0.D0
 382       wturn6=0.D0
 383       wang=0.D0
 384       wtor=0.D0
 385       wtor_d=0.D0
 386
 387 c     Prepare to freeze backbone
 388       do i=1,nres
 389         mask_phi(i)=0
 390         mask_theta(i)=0
 391         mask_side(i)=1
 392       enddo
 393
 394 c     Minimize the side-chains
 395       mask_r=.true.
 396       call geom_to_var(nvar,var)
 397       call minimize(etot,var,iretcode,nfun)
 398       call var_to_geom(nvar,var)
 399       mask_r=.false.
 400
 401 c     Put the original weights back and calculate the full energy
 402       wscp=orig_w(2)
 403       welec=orig_w(3)
 404       wcorr=orig_w(4)
 405       wcorr5=orig_w(5)
 406       wcorr6=orig_w(6)
 407       wel_loc=orig_w(7)
 408       wturn3=orig_w(8)
 409       wturn4=orig_w(9)
 410       wturn6=orig_w(10)
 411       wang=orig_w(11)
 412       wtor=orig_w(13)
 413       wtor_d=orig_w(14)
 414
 415       call chainbuild
 416       call etotal(energy)
 417       etot=energy(0)
 418
 419       return
 420       end
 421
 422 c-------------------------------------------------------------
 423       subroutine minimize_sc1(etot,x,iretcode,nfun)
 424       implicit real*8 (a-h,o-z)
 425       include 'DIMENSIONS'
 426       parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
 427       include 'COMMON.IOUNITS'
 428       include 'COMMON.VAR'
 429       include 'COMMON.GEO'
 430       include 'COMMON.MINIM'
 431       common /srutu/ icall
 432       dimension iv(liv)
 433       double precision minval,x(maxvar),d(maxvar),v(1:lv),xx(maxvar)
 434       double precision energia(0:n_ene)
 435       external func,gradient,fdum
 436       external func_restr1,grad_restr1
 437       logical not_done,change,reduce
 438       common /przechowalnia/ v
 439
 440       call deflt(2,iv,liv,lv,v)
 441 * 12 means fresh start, dont call deflt
 442       iv(1)=12
 443 * max num of fun calls
 444       if (maxfun.eq.0) maxfun=500
 445       iv(17)=maxfun
 446 * max num of iterations
 447       if (maxmin.eq.0) maxmin=1000
 448       iv(18)=maxmin
 449 * controls output
 450       iv(19)=2
 451 * selects output unit
 452 c     iv(21)=iout
 453       iv(21)=0
 454 * 1 means to print out result
 455       iv(22)=0
 456 * 1 means to print out summary stats
 457       iv(23)=0
 458 * 1 means to print initial x and d
 459       iv(24)=0
 460 * min val for v(radfac) default is 0.1
 461       v(24)=0.1D0
 462 * max val for v(radfac) default is 4.0
 463       v(25)=2.0D0
 464 c     v(25)=4.0D0
 465 * check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
 466 * the sumsl default is 0.1
 467       v(26)=0.1D0
 468 * false conv if (act fnctn decrease) .lt. v(34)
 469 * the sumsl default is 100*machep
 470       v(34)=v(34)/100.0D0
 471 * absolute convergence
 472       if (tolf.eq.0.0D0) tolf=1.0D-4
 473       v(31)=tolf
 474 * relative convergence
 475       if (rtolf.eq.0.0D0) rtolf=1.0D-4
 476       v(32)=rtolf
 477 * controls initial step size
 478        v(35)=1.0D-1
 479 * large vals of d correspond to small components of step
 480       do i=1,nphi
 481         d(i)=1.0D-1
 482       enddo
 483       do i=nphi+1,nvar
 484         d(i)=1.0D-1
 485       enddo
 486       IF (mask_r) THEN
 487        call x2xx(x,xx,nvar_restr)
 488        call sumsl(nvar_restr,d,xx,func_restr1,grad_restr1,
 489      &                    iv,liv,lv,v,idum,rdum,fdum)
 490        call xx2x(x,xx)
 491       ELSE
 492        call sumsl(nvar,d,x,func,gradient,iv,liv,lv,v,idum,rdum,fdum)
 493       ENDIF
 494       etot=v(10)
 495       iretcode=iv(1)
 496       nfun=iv(6)
 497
 498       return
 499       end
 500 ************************************************************************
 501       subroutine func_restr1(n,x,nf,f,uiparm,urparm,ufparm)
 502       implicit real*8 (a-h,o-z)
 503       include 'DIMENSIONS'
 504       include 'COMMON.DERIV'
 505       include 'COMMON.IOUNITS'
 506       include 'COMMON.GEO'
 507       include 'COMMON.FFIELD'
 508       include 'COMMON.INTERACT'
 509       include 'COMMON.TIME1'
 510       common /chuju/ jjj
 511       double precision energia(0:n_ene),evdw,escloc
 512       integer jjj
 513       double precision ufparm,e1,e2
 514       external ufparm
 515       integer uiparm(1)
 516       real*8 urparm(1)
 517       dimension x(maxvar)
 518       nfl=nf
 519       icg=mod(nf,2)+1
 520
 521 #ifdef OSF
 522 c     Intercept NaNs in the coordinates, before calling etotal
 523       x_sum=0.D0
 524       do i=1,n
 525         x_sum=x_sum+x(i)
 526       enddo
 527       FOUND_NAN=.false.
 528       if (x_sum.ne.x_sum) then
 529         write(iout,*)"   *** func_restr1 : Found NaN in coordinates"
 530         f=1.0D+73
 531         FOUND_NAN=.true.
 532         return
 533       endif
 534 #endif
 535
 536       call var_to_geom_restr(n,x)
 537       call zerograd
 538       call chainbuild
 539 cd    write (iout,*) 'ETOTAL called from FUNC'
 540       call egb1(evdw)
 541       call esc(escloc)
 542       f=wsc*evdw+wscloc*escloc
 543 cd      call etotal(energia(0))
 544 cd      f=wsc*energia(1)+wscloc*energia(12)
 545 cd      print *,f,evdw,escloc,energia(0)
 546 C
 547 C Sum up the components of the Cartesian gradient.
 548 C
 549       do i=1,nct
 550         do j=1,3
 551           gradx(j,i,icg)=wsc*gvdwx(j,i)
 552         enddo
 553       enddo
 554
 555       return
 556       end
 557 c-------------------------------------------------------
 558       subroutine grad_restr1(n,x,nf,g,uiparm,urparm,ufparm)
 559       implicit real*8 (a-h,o-z)
 560       include 'DIMENSIONS'
 561       include 'COMMON.CHAIN'
 562       include 'COMMON.DERIV'
 563       include 'COMMON.VAR'
 564       include 'COMMON.INTERACT'
 565       include 'COMMON.FFIELD'
 566       include 'COMMON.IOUNITS'
 567       external ufparm
 568       integer uiparm(1)
 569       double precision urparm(1)
 570       dimension x(maxvar),g(maxvar)
 571
 572       icg=mod(nf,2)+1
 573       if (nf-nfl+1) 20,30,40
 574    20 call func_restr1(n,x,nf,f,uiparm,urparm,ufparm)
 575 c     write (iout,*) 'grad 20'
 576       if (nf.eq.0) return
 577       goto 40
 578    30 call var_to_geom_restr(n,x)
 579       call chainbuild
 580 C
 581 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
 582 C
 583    40 call cartder
 584 C
 585 C Convert the Cartesian gradient into internal-coordinate gradient.
 586 C
 587
 588       ig=0
 589       ind=nres-2
 590       do i=2,nres-2
 591        IF (mask_phi(i+2).eq.1) THEN
 592         gphii=0.0D0
 593         do j=i+1,nres-1
 594           ind=ind+1
 595           do k=1,3
 596             gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
 597             gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)
 598           enddo
 599         enddo
 600         ig=ig+1
 601         g(ig)=gphii
 602        ELSE
 603         ind=ind+nres-1-i
 604        ENDIF
 605       enddo
 606
 607
 608       ind=0
 609       do i=1,nres-2
 610        IF (mask_theta(i+2).eq.1) THEN
 611         ig=ig+1
 612         gthetai=0.0D0
 613         do j=i+1,nres-1
 614           ind=ind+1
 615           do k=1,3
 616             gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
 617             gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
 618           enddo
 619         enddo
 620         g(ig)=gthetai
 621        ELSE
 622         ind=ind+nres-1-i
 623        ENDIF
 624       enddo
 625
 626       do i=2,nres-1
 627         if (itype(i).ne.10) then
 628          IF (mask_side(i).eq.1) THEN
 629           ig=ig+1
 630           galphai=0.0D0
 631           do k=1,3
 632             galphai=galphai+dxds(k,i)*gradx(k,i,icg)
 633           enddo
 634           g(ig)=galphai
 635          ENDIF
 636         endif
 637       enddo
 638
 639
 640       do i=2,nres-1
 641         if (itype(i).ne.10) then
 642          IF (mask_side(i).eq.1) THEN
 643           ig=ig+1
 644           gomegai=0.0D0
 645           do k=1,3
 646             gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
 647           enddo
 648           g(ig)=gomegai
 649          ENDIF
 650         endif
 651       enddo
 652
 653 C
 654 C Add the components corresponding to local energy terms.
 655 C
 656
 657       ig=0
 658       igall=0
 659       do i=4,nres
 660         igall=igall+1
 661         if (mask_phi(i).eq.1) then
 662           ig=ig+1
 663           g(ig)=g(ig)+gloc(igall,icg)
 664         endif
 665       enddo
 666
 667       do i=3,nres
 668         igall=igall+1
 669         if (mask_theta(i).eq.1) then
 670           ig=ig+1
 671           g(ig)=g(ig)+gloc(igall,icg)
 672         endif
 673       enddo
 674
 675       do ij=1,2
 676       do i=2,nres-1
 677         if (itype(i).ne.10) then
 678           igall=igall+1
 679           if (mask_side(i).eq.1) then
 680             ig=ig+1
 681             g(ig)=g(ig)+gloc(igall,icg)
 682           endif
 683         endif
 684       enddo
 685       enddo
 686
 687 cd      do i=1,ig
 688 cd        write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
 689 cd      enddo
 690       return
 691       end
 692 C-----------------------------------------------------------------------------
 693       subroutine egb1(evdw)
 694 C
 695 C This subroutine calculates the interaction energy of nonbonded side chains
 696 C assuming the Gay-Berne potential of interaction.
 697 C
 698       implicit real*8 (a-h,o-z)
 699       include 'DIMENSIONS'
 700       include 'COMMON.GEO'
 701       include 'COMMON.VAR'
 702       include 'COMMON.LOCAL'
 703       include 'COMMON.CHAIN'
 704       include 'COMMON.DERIV'
 705       include 'COMMON.NAMES'
 706       include 'COMMON.INTERACT'
 707       include 'COMMON.IOUNITS'
 708       include 'COMMON.CALC'
 709       include 'COMMON.CONTROL'
 710       logical lprn
 711       evdw=0.0D0
 712 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
 713       evdw=0.0D0
 714       lprn=.false.
 715 c     if (icall.eq.0) lprn=.true.
 716       ind=0
 717       do i=iatsc_s,iatsc_e
 718
 719
 720         itypi=itype(i)
 721         itypi1=itype(i+1)
 722         xi=c(1,nres+i)
 723         yi=c(2,nres+i)
 724         zi=c(3,nres+i)
 725         dxi=dc_norm(1,nres+i)
 726         dyi=dc_norm(2,nres+i)
 727         dzi=dc_norm(3,nres+i)
 728         dsci_inv=dsc_inv(itypi)
 729 C
 730 C Calculate SC interaction energy.
 731 C
 732         do iint=1,nint_gr(i)
 733           do j=istart(i,iint),iend(i,iint)
 734           IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
 735             ind=ind+1
 736             itypj=itype(j)
 737             dscj_inv=dsc_inv(itypj)
 738             sig0ij=sigma(itypi,itypj)
 739             chi1=chi(itypi,itypj)
 740             chi2=chi(itypj,itypi)
 741             chi12=chi1*chi2
 742             chip1=chip(itypi)
 743             chip2=chip(itypj)
 744             chip12=chip1*chip2
 745             alf1=alp(itypi)
 746             alf2=alp(itypj)
 747             alf12=0.5D0*(alf1+alf2)
 748 C For diagnostics only!!!
 749 c           chi1=0.0D0
 750 c           chi2=0.0D0
 751 c           chi12=0.0D0
 752 c           chip1=0.0D0
 753 c           chip2=0.0D0
 754 c           chip12=0.0D0
 755 c           alf1=0.0D0
 756 c           alf2=0.0D0
 757 c           alf12=0.0D0
 758             xj=c(1,nres+j)-xi
 759             yj=c(2,nres+j)-yi
 760             zj=c(3,nres+j)-zi
 761             dxj=dc_norm(1,nres+j)
 762             dyj=dc_norm(2,nres+j)
 763             dzj=dc_norm(3,nres+j)
 764             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
 765             rij=dsqrt(rrij)
 766 C Calculate angle-dependent terms of energy and contributions to their
 767 C derivatives.
 768             call sc_angular
 769             sigsq=1.0D0/sigsq
 770             sig=sig0ij*dsqrt(sigsq)
 771             rij_shift=1.0D0/rij-sig+sig0ij
 772 C I hate to put IF's in the loops, but here don't have another choice!!!!
 773             if (rij_shift.le.0.0D0) then
 774               evdw=1.0D20
 775 cd              write (iout,'(2(a3,i3,2x),17(0pf7.3))')
 776 cd     &        restyp(itypi),i,restyp(itypj),j,
 777 cd     &        rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
 778               return
 779             endif
 780             sigder=-sig*sigsq
 781 c---------------------------------------------------------------
 782             rij_shift=1.0D0/rij_shift
 783             fac=rij_shift**expon
 784             e1=fac*fac*aa(itypi,itypj)
 785             e2=fac*bb(itypi,itypj)
 786             evdwij=eps1*eps2rt*eps3rt*(e1+e2)
 787             eps2der=evdwij*eps3rt
 788             eps3der=evdwij*eps2rt
 789             evdwij=evdwij*eps2rt*eps3rt
 790             evdw=evdw+evdwij
 791             if (lprn) then
 792             sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
 793             epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
 794 cd            write (iout,'(2(a3,i3,2x),17(0pf7.3))')
 795 cd     &        restyp(itypi),i,restyp(itypj),j,
 796 cd     &        epsi,sigm,chi1,chi2,chip1,chip2,
 797 cd     &        eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
 798 cd     &        om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
 799 cd     &        evdwij
 800             endif
 801
 802             if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
 803      &                        'evdw',i,j,evdwij
 804
 805 C Calculate gradient components.
 806             e1=e1*eps1*eps2rt**2*eps3rt**2
 807             fac=-expon*(e1+evdwij)*rij_shift
 808             sigder=fac*sigder
 809             fac=rij*fac
 810 C Calculate the radial part of the gradient
 811             gg(1)=xj*fac
 812             gg(2)=yj*fac
 813             gg(3)=zj*fac
 814 C Calculate angular part of the gradient.
 815             call sc_grad
 816           ENDIF
 817           enddo      ! j
 818         enddo        ! iint
 819       enddo          ! i
 820       end
 821 C-----------------------------------------------------------------------------