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