source/unres/src-HCD-5D/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       sideonly=.true.
  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       sideonly=.false.
 156       mask_side=1
 157 crc      n_fun=n_fun+1
 158 ct      write (iout,*) 'sc_local time= ',MPI_WTIME()-wtime
 159       call chainbuild_extconf
 160       return
 161       end
 162
 163 c-------------------------------------------------------------
 164
 165       subroutine single_sc_move(res_pick,n_maxtry,e_drop,
 166      +     n_steps,n_fun,e_sc)
 167 c     Perturb one side-chain (res_pick) and minimize the
 168 c     neighbouring region, keeping all CA's and non-neighbouring
 169 c     side-chains fixed
 170 c     Try until e_drop energy improvement is achieved, or n_maxtry
 171 c     attempts have been made
 172 c     At the start, e_sc should contain the side-chain-only energy(0)
 173 c     nsteps and nfun for this move are ADDED to n_steps and n_fun
 174 crc      implicit none
 175
 176 c     Includes
 177       implicit real*8 (a-h,o-z)
 178       include 'DIMENSIONS'
 179       include 'COMMON.VAR'
 180       include 'COMMON.INTERACT'
 181       include 'COMMON.CHAIN'
 182       include 'COMMON.MINIM'
 183       include 'COMMON.FFIELD'
 184       include 'COMMON.IOUNITS'
 185
 186 c     External functions
 187       double precision dist
 188       external dist
 189
 190 c     Input arguments
 191       integer res_pick,n_maxtry
 192       double precision e_drop
 193
 194 c     Input/Output arguments
 195       integer n_steps,n_fun
 196       double precision e_sc
 197
 198 c     Local variables
 199       logical fail
 200       integer i,j
 201       integer nres_moved
 202       integer iretcode,loc_nfun,orig_maxfun,n_try
 203       double precision sc_dist,sc_dist_cutoff
 204       double precision energy(0:n_ene),orig_e,cur_e
 205       double precision evdw,escloc
 206       double precision cur_alph(2:nres-1),cur_omeg(2:nres-1)
 207       double precision var(maxvar)
 208
 209       double precision orig_theta(1:nres),orig_phi(1:nres),
 210      +     orig_alph(1:nres),orig_omeg(1:nres)
 211
 212
 213 c     Define what is meant by "neighbouring side-chain"
 214       sc_dist_cutoff=5.0D0
 215
 216 c     Don't do glycine or ends
 217       i=itype(res_pick)
 218       if (i.eq.10 .or. i.eq.ntyp1) return
 219
 220 c     Freeze everything (later will relax only selected side-chains)
 221       mask_r=.true.
 222       do i=1,nres
 223         mask_phi(i)=0
 224         mask_theta(i)=0
 225         mask_side(i)=0
 226       enddo
 227
 228 c     Find the neighbours of the side-chain to move
 229 c     and save initial variables
 230 crc      orig_e=e_sc
 231 crc      cur_e=orig_e
 232       nres_moved=0
 233       do i=2,nres-1
 234 c     Don't do glycine (itype(j)==10)
 235         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 236           sc_dist=dist(nres+i,nres+res_pick)
 237         else
 238           sc_dist=sc_dist_cutoff
 239         endif
 240         if (sc_dist.lt.sc_dist_cutoff) then
 241           nres_moved=nres_moved+1
 242           mask_side(i)=1
 243           cur_alph(i)=alph(i)
 244           cur_omeg(i)=omeg(i)
 245         endif
 246       enddo
 247
 248       call chainbuild_extconf
 249       call egb1(evdw)
 250       call esc(escloc)
 251       e_sc=wsc*evdw+wscloc*escloc
 252 c      write (iout,*) "sc_move: e_sc",e_sc
 253 cd      call etotal(energy)
 254 cd      print *,'new       ',(energy(k),k=0,n_ene)
 255       orig_e=e_sc
 256       cur_e=orig_e
 257
 258       n_try=0
 259       do while (n_try.lt.n_maxtry .and. orig_e-cur_e.lt.e_drop)
 260 c     Move the selected residue (don't worry if it fails)
 261         call gen_side(iabs(itype(res_pick)),theta(res_pick+1),
 262      +       alph(res_pick),omeg(res_pick),fail)
 263
 264 c     Minimize the side-chains starting from the new arrangement
 265         call geom_to_var(nvar,var)
 266         orig_maxfun=maxfun
 267         maxfun=7
 268
 269 crc        do i=1,nres
 270 crc          orig_theta(i)=theta(i)
 271 crc          orig_phi(i)=phi(i)
 272 crc          orig_alph(i)=alph(i)
 273 crc          orig_omeg(i)=omeg(i)
 274 crc        enddo
 275
 276         call minimize_sc1(e_sc,var,iretcode,loc_nfun)
 277 c        write (iout,*) "n_try",n_try
 278 c        write (iout,*) "sc_move after minimze_sc1 e_sc",e_sc
 279 cv        write(*,'(2i3,2f12.5,2i3)')
 280 cv     &       res_pick,nres_moved,orig_e,e_sc-cur_e,
 281 cv     &       iretcode,loc_nfun
 282
 283 c$$$        if (iretcode.eq.8) then
 284 c$$$          write(iout,*)'Coordinates just after code 8'
 285 c$$$          call chainbuild
 286 c$$$          call all_varout
 287 c$$$          call flush(iout)
 288 c$$$          do i=1,nres
 289 c$$$            theta(i)=orig_theta(i)
 290 c$$$            phi(i)=orig_phi(i)
 291 c$$$            alph(i)=orig_alph(i)
 292 c$$$            omeg(i)=orig_omeg(i)
 293 c$$$          enddo
 294 c$$$          write(iout,*)'Coordinates just before code 8'
 295 c$$$          call chainbuild
 296 c$$$          call all_varout
 297 c$$$          call flush(iout)
 298 c$$$        endif
 299
 300         n_fun=n_fun+loc_nfun
 301         maxfun=orig_maxfun
 302         call var_to_geom(nvar,var)
 303
 304 c     If a lower energy was found, update the current structure...
 305         if (e_sc.lt.cur_e) then
 306 cv              call chainbuild
 307 cv              call etotal(energy)
 308 cd              call egb1(evdw)
 309 cd              call esc(escloc)
 310 cd              e_sc1=wsc*evdw+wscloc*escloc
 311 cd              print *,'     new',e_sc1,energy(0)
 312 cv              print *,'new       ',energy(0)
 313 cd              call enerprint(energy(0))
 314           cur_e=e_sc
 315           do i=2,nres-1
 316             if (mask_side(i).eq.1) then
 317               cur_alph(i)=alph(i)
 318               cur_omeg(i)=omeg(i)
 319             endif
 320           enddo
 321         else
 322 c     ...else revert to the previous one
 323           e_sc=cur_e
 324           do i=2,nres-1
 325             if (mask_side(i).eq.1) then
 326               alph(i)=cur_alph(i)
 327               omeg(i)=cur_omeg(i)
 328             endif
 329           enddo
 330         endif
 331         n_try=n_try+1
 332
 333       enddo
 334       n_steps=n_steps+n_try
 335
 336 c     Reset the minimization mask_r to false
 337       mask_r=.false.
 338
 339       return
 340       end
 341 c-------------------------------------------------------------
 342       subroutine minimize_sc1(etot,x,iretcode,nfun)
 343 #ifdef LBFGS_SC
 344       use minima
 345       use inform
 346       use output
 347       use iounit
 348       use scales
 349 #endif
 350       implicit real*8 (a-h,o-z)
 351       include 'DIMENSIONS'
 352 #ifndef LBFGS_SC
 353 c      parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
 354       parameter(max_sc_move=10)
 355       parameter (liv=60,lv=(77+2*max_sc_move*(2*max_sc_move+17)/2))
 356 #endif
 357       include 'COMMON.IOUNITS'
 358       include 'COMMON.VAR'
 359       include 'COMMON.GEO'
 360       include 'COMMON.MINIM'
 361       common /srutu/ icall
 362       double precision x(maxvar),d(maxvar),xx(maxvar)
 363       double precision energia(0:n_ene)
 364 #ifdef LBFGS_SC
 365       integer nvar_restr
 366       common /zmienne/ nvar_restr
 367       double precision grdmin
 368       double precision funcgrad_restr1
 369       external funcgrad_restr1
 370       external optsave
 371 #else
 372       external func,gradient,fdum
 373       external func_restr1,grad_restr1
 374       logical not_done,change,reduce
 375       dimension iv(liv)
 376       double precision v(1:lv)
 377       common /przechowalnia/ v
 378 #endif
 379 #ifdef LBFGS_SC
 380       maxiter=7
 381       coordtype='RIGIDBODY'
 382       grdmin=tolf
 383       jout=iout
 384 c      jprint=print_min_stat
 385       jprint=0
 386       iwrite=0
 387       if (.not. allocated(scale))  allocate (scale(nvar))
 388 c
 389 c     set scaling parameter for function and derivative values;
 390 c     use square root of median eigenvalue of typical Hessian
 391 c
 392       call x2xx(x,xx,nvar_restr)
 393       set_scale = .true.
 394 c      nvar = 0
 395       do i = 1, nvar_restr
 396 c         if (use(i)) then
 397 c            do j = 1, 3
 398 c               nvar = nvar + 1
 399                scale(i) = 12.0d0
 400 c            end do
 401 c         end if
 402       end do
 403 c      write (iout,*) "Calling lbfgs"
 404       call lbfgs (nvar_restr,xx,etot,grdmin,funcgrad_restr1,optsave)
 405       deallocate(scale)
 406 c      write (iout,*) "After lbfgs"
 407       call xx2x(x,xx)
 408 #else
 409       call deflt(2,iv,liv,lv,v)
 410 * 12 means fresh start, dont call deflt
 411       iv(1)=12
 412 * max num of fun calls
 413       if (maxfun.eq.0) maxfun=500
 414       iv(17)=maxfun
 415 * max num of iterations
 416       if (maxmin.eq.0) maxmin=1000
 417       iv(18)=maxmin
 418 * controls output
 419       iv(19)=2
 420 * selects output unit
 421       iv(21)=0
 422 c      iv(21)=0
 423 * 1 means to print out result
 424       iv(22)=0
 425 * 1 means to print out summary stats
 426       iv(23)=0
 427 * 1 means to print initial x and d
 428       iv(24)=0
 429 * min val for v(radfac) default is 0.1
 430       v(24)=0.1D0
 431 * max val for v(radfac) default is 4.0
 432       v(25)=2.0D0
 433 c     v(25)=4.0D0
 434 * check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
 435 * the sumsl default is 0.1
 436       v(26)=0.1D0
 437 * false conv if (act fnctn decrease) .lt. v(34)
 438 * the sumsl default is 100*machep
 439       v(34)=v(34)/100.0D0
 440 * absolute convergence
 441       if (tolf.eq.0.0D0) tolf=1.0D-4
 442       v(31)=tolf
 443 * relative convergence
 444       if (rtolf.eq.0.0D0) rtolf=1.0D-4
 445       v(32)=rtolf
 446 * controls initial step size
 447        v(35)=1.0D-1
 448 * large vals of d correspond to small components of step
 449       do i=1,nphi
 450         d(i)=1.0D-1
 451       enddo
 452       do i=nphi+1,nvar
 453         d(i)=1.0D-1
 454       enddo
 455       IF (mask_r) THEN
 456        call x2xx(x,xx,nvar_restr)
 457        call sumsl(nvar_restr,d,xx,func_restr1,grad_restr1,
 458      &                    iv,liv,lv,v,idum,rdum,fdum)
 459        call xx2x(x,xx)
 460       ELSE
 461 c       call sumsl(nvar,d,x,func,gradient,iv,liv,lv,v,idum,rdum,fdum)
 462       ENDIF
 463       etot=v(10)
 464       iretcode=iv(1)
 465       nfun=iv(6)
 466 #endif
 467       return
 468       end
 469 #ifdef LBFGS_SC
 470       double precision function funcgrad_restr1(x,g)
 471       implicit real*8 (a-h,o-z)
 472       include 'DIMENSIONS'
 473       include 'COMMON.DERIV'
 474       include 'COMMON.IOUNITS'
 475       include 'COMMON.GEO'
 476       include 'COMMON.FFIELD'
 477       include 'COMMON.INTERACT'
 478       include 'COMMON.TIME1'
 479       include 'COMMON.CHAIN'
 480       include 'COMMON.VAR'
 481       integer nvar_restr
 482       common /zmienne/ nvar_restr
 483       double precision energia(0:n_ene),evdw,escloc
 484       double precision ufparm,e1,e2
 485       dimension x(maxvar),g(maxvar),gg(maxvar)
 486 #ifdef OSF
 487 c     Intercept NaNs in the coordinates, before calling etotal
 488       x_sum=0.D0
 489       do i=1,nvar_restr
 490         x_sum=x_sum+x(i)
 491       enddo
 492       FOUND_NAN=.false.
 493       if (x_sum.ne.x_sum) then
 494         write(iout,*)"   *** func_restr1 : Found NaN in coordinates"
 495         f=1.0D+73
 496         FOUND_NAN=.true.
 497         return
 498       endif
 499 #else
 500       FOUND_NAN=.false.
 501       do i=1,nvar_restr
 502       if (isnan(x(i))) then
 503         FOUND_NAN=.true.
 504         f=1.0D+73
 505         funcgrad_restr1=f
 506         write (iout,*) "NaN in coordinates"
 507         return
 508       endif
 509       enddo
 510 #endif
 511
 512 c      write (iout,*) "nvar_restr",nvar_restr
 513 c      write (iout,*) "x",(x(i),i=1,nvar_restr)
 514       call var_to_geom_restr(nvar_restr,x)
 515       call zerograd
 516       call chainbuild_extconf
 517 cd    write (iout,*) 'ETOTAL called from FUNC'
 518       call egb1(evdw)
 519       call esc(escloc)
 520       f=wsc*evdw+wscloc*escloc
 521 c      write (iout,*) "evdw",evdw," escloc",escloc
 522       if (isnan(f)) then
 523         f=1.0d20
 524         funcgrad_restr1=f
 525         return
 526       endif
 527       funcgrad_restr1=f
 528 c      write (iout,*) "f",f
 529 cd      call etotal(energia(0))
 530 cd      f=wsc*energia(1)+wscloc*energia(12)
 531 cd      print *,f,evdw,escloc,energia(0)
 532 C
 533 C Sum up the components of the Cartesian gradient.
 534 C
 535       do i=1,nct
 536         do j=1,3
 537           gradx(j,i,icg)=wsc*gvdwx(j,i)+wscloc*gsclocx(j,i)
 538         enddo
 539       enddo
 540
 541 C
 542 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
 543 C
 544       call cart2intgrad(nvar,gg)
 545 C
 546 C Convert the Cartesian gradient into internal-coordinate gradient.
 547 C
 548
 549       ig=0
 550       do i=2,nres-1
 551         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 552          IF (mask_side(i).eq.1) THEN
 553           ig=ig+1
 554           g(ig)=gg(ialph(i,1))
 555 c          write (iout,*) "i",i," ig",ig," ialph",ialph(i,1)
 556 c          write (iout,*) "g",g(ig)," gg",gg(ialph(i,1))
 557          ENDIF
 558         endif
 559       enddo
 560       do i=2,nres-1
 561         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 562          IF (mask_side(i).eq.1) THEN
 563           ig=ig+1
 564           g(ig)=gg(ialph(i,1)+nside)
 565 c          write (iout,*) "i",i," ig",ig," ialph",ialph(i,1)+nside
 566 c          write (iout,*) "g",g(ig)," gg",gg(ialph(i,1)+nside)
 567          ENDIF
 568         endif
 569       enddo
 570
 571 C
 572 C Add the components corresponding to local energy terms.
 573 C
 574
 575       ig=0
 576       igall=0
 577       do i=4,nres
 578         igall=igall+1
 579         if (mask_phi(i).eq.1) then
 580           ig=ig+1
 581           g(ig)=g(ig)+gloc(igall,icg)
 582         endif
 583       enddo
 584
 585       do i=3,nres
 586         igall=igall+1
 587         if (mask_theta(i).eq.1) then
 588           ig=ig+1
 589           g(ig)=g(ig)+gloc(igall,icg)
 590         endif
 591       enddo
 592
 593       do ij=1,2
 594       do i=2,nres-1
 595         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 596           igall=igall+1
 597           if (mask_side(i).eq.1) then
 598             ig=ig+1
 599             g(ig)=g(ig)+gloc(igall,icg)
 600 c            write (iout,*) "ij",ij," i",i," ig",ig," igall",igall
 601 c            write (iout,*) "gloc",gloc(igall,icg)," g",g(ig)
 602           endif
 603         endif
 604       enddo
 605       enddo
 606
 607 cd      do i=1,ig
 608 cd        write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
 609 cd      enddo
 610       return
 611       end
 612 #else
 613 ************************************************************************
 614       subroutine func_restr1(n,x,nf,f,uiparm,urparm,ufparm)
 615       implicit real*8 (a-h,o-z)
 616       include 'DIMENSIONS'
 617       include 'COMMON.DERIV'
 618       include 'COMMON.IOUNITS'
 619       include 'COMMON.GEO'
 620       include 'COMMON.FFIELD'
 621       include 'COMMON.INTERACT'
 622       include 'COMMON.TIME1'
 623       double precision energia(0:n_ene),evdw,escloc
 624       double precision ufparm,e1,e2
 625       external ufparm
 626       integer uiparm(1)
 627       real*8 urparm(1)
 628       dimension x(maxvar)
 629       nfl=nf
 630       icg=mod(nf,2)+1
 631
 632 #ifdef OSF
 633 c     Intercept NaNs in the coordinates, before calling etotal
 634       x_sum=0.D0
 635       do i=1,n
 636         x_sum=x_sum+x(i)
 637       enddo
 638       FOUND_NAN=.false.
 639       if (x_sum.ne.x_sum) then
 640         write(iout,*)"   *** func_restr1 : Found NaN in coordinates"
 641         f=1.0D+73
 642         FOUND_NAN=.true.
 643         return
 644       endif
 645 #endif
 646
 647       call var_to_geom_restr(n,x)
 648       call zerograd
 649       call chainbuild_extconf
 650 cd    write (iout,*) 'ETOTAL called from FUNC'
 651       call egb1(evdw)
 652       call esc(escloc)
 653       f=wsc*evdw+wscloc*escloc
 654 c      write (iout,*) "f",f
 655 cd      call etotal(energia(0))
 656 cd      f=wsc*energia(1)+wscloc*energia(12)
 657 cd      print *,f,evdw,escloc,energia(0)
 658 C
 659 C Sum up the components of the Cartesian gradient.
 660 C
 661       do i=1,nct
 662         do j=1,3
 663           gradx(j,i,icg)=wsc*gvdwx(j,i)+wscloc*gsclocx(j,i)
 664         enddo
 665       enddo
 666
 667       return
 668       end
 669 c-------------------------------------------------------
 670       subroutine grad_restr1(n,x,nf,g,uiparm,urparm,ufparm)
 671       implicit real*8 (a-h,o-z)
 672       include 'DIMENSIONS'
 673       include 'COMMON.CHAIN'
 674       include 'COMMON.DERIV'
 675       include 'COMMON.VAR'
 676       include 'COMMON.INTERACT'
 677       include 'COMMON.FFIELD'
 678       include 'COMMON.IOUNITS'
 679       external ufparm
 680       integer uiparm(1)
 681       double precision urparm(1)
 682       dimension x(maxvar),g(maxvar),gg(maxvar)
 683
 684       icg=mod(nf,2)+1
 685       if (nf-nfl+1) 20,30,40
 686    20 call func_restr1(n,x,nf,f,uiparm,urparm,ufparm)
 687 c     write (iout,*) 'grad 20'
 688       if (nf.eq.0) return
 689       goto 40
 690    30 call var_to_geom_restr(n,x)
 691       call chainbuild_extconf
 692 C
 693 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
 694 C
 695    40 call cart2intgrad(nvar,gg)
 696 C
 697 C Convert the Cartesian gradient into internal-coordinate gradient.
 698 C
 699
 700       ig=0
 701       ind=nres-2
 702       do i=2,nres-2
 703        IF (mask_phi(i+2).eq.1) THEN
 704         ig=ig+1
 705         g(ig)=gg(i-1)
 706        ENDIF
 707       enddo
 708
 709
 710       do i=1,nres-2
 711        IF (mask_theta(i+2).eq.1) THEN
 712         ig=ig+1
 713         g(ig)=gg(nphi+i)
 714        ENDIF
 715       enddo
 716
 717       do i=2,nres-1
 718         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 719          IF (mask_side(i).eq.1) THEN
 720           ig=ig+1
 721           g(ig)=gg(ialph(i,1))
 722 c          write (iout,*) "i",i," ig",ig," ialph",ialph(i,1)
 723 c          write (iout,*) "g",g(ig)," gg",gg(ialph(i,1))
 724          ENDIF
 725         endif
 726       enddo
 727
 728
 729       do i=2,nres-1
 730         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 731          IF (mask_side(i).eq.1) THEN
 732           ig=ig+1
 733           g(ig)=gg(ialph(i,1)+nside)
 734 c          write (iout,*) "i",i," ig",ig," ialph",ialph(i,1)+nside
 735 c          write (iout,*) "g",g(ig)," gg",gg(ialph(i,1)+nside)
 736          ENDIF
 737         endif
 738       enddo
 739
 740 C
 741 C Add the components corresponding to local energy terms.
 742 C
 743
 744       ig=0
 745       igall=0
 746       do i=4,nres
 747         igall=igall+1
 748         if (mask_phi(i).eq.1) then
 749           ig=ig+1
 750           g(ig)=g(ig)+gloc(igall,icg)
 751         endif
 752       enddo
 753
 754       do i=3,nres
 755         igall=igall+1
 756         if (mask_theta(i).eq.1) then
 757           ig=ig+1
 758           g(ig)=g(ig)+gloc(igall,icg)
 759         endif
 760       enddo
 761
 762       do ij=1,2
 763       do i=2,nres-1
 764         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 765           igall=igall+1
 766           if (mask_side(i).eq.1) then
 767             ig=ig+1
 768             g(ig)=g(ig)+gloc(igall,icg)
 769 c            write (iout,*) "ij",ij," i",i," ig",ig," igall",igall
 770 c            write (iout,*) "gloc",gloc(igall,icg)," g",g(ig)
 771           endif
 772         endif
 773       enddo
 774       enddo
 775
 776 cd      do i=1,ig
 777 cd        write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
 778 cd      enddo
 779       return
 780       end
 781 #endif
 782 C-----------------------------------------------------------------------------
 783       subroutine egb1(evdw)
 784 C
 785 C This subroutine calculates the interaction energy of nonbonded side chains
 786 C assuming the Gay-Berne potential of interaction.
 787 C
 788       implicit real*8 (a-h,o-z)
 789       include 'DIMENSIONS'
 790       include 'COMMON.GEO'
 791       include 'COMMON.VAR'
 792       include 'COMMON.LOCAL'
 793       include 'COMMON.CHAIN'
 794       include 'COMMON.DERIV'
 795       include 'COMMON.NAMES'
 796       include 'COMMON.INTERACT'
 797       include 'COMMON.IOUNITS'
 798       include 'COMMON.CALC'
 799       include 'COMMON.CONTROL'
 800       logical lprn
 801       evdw=0.0D0
 802 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
 803       evdw=0.0D0
 804       lprn=.false.
 805 c     if (icall.eq.0) lprn=.true.
 806       ind=0
 807 c      do i=iatsc_s,iatsc_e
 808       do i=nnt,nct
 809
 810
 811         itypi=iabs(itype(i))
 812         if (itypi.eq.ntyp1 .or. mask_side(i).eq.0) cycle
 813         itypi1=iabs(itype(i+1))
 814         xi=c(1,nres+i)
 815         yi=c(2,nres+i)
 816         zi=c(3,nres+i)
 817           xi=mod(xi,boxxsize)
 818           if (xi.lt.0) xi=xi+boxxsize
 819           yi=mod(yi,boxysize)
 820           if (yi.lt.0) yi=yi+boxysize
 821           zi=mod(zi,boxzsize)
 822           if (zi.lt.0) zi=zi+boxzsize
 823        if ((zi.gt.bordlipbot)
 824      &.and.(zi.lt.bordliptop)) then
 825 C the energy transfer exist
 826         if (zi.lt.buflipbot) then
 827 C what fraction I am in
 828          fracinbuf=1.0d0-
 829      &        ((positi-bordlipbot)/lipbufthick)
 830 C lipbufthick is thickenes of lipid buffore
 831          sslipi=sscalelip(fracinbuf)
 832          ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
 833         elseif (zi.gt.bufliptop) then
 834          fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
 835          sslipi=sscalelip(fracinbuf)
 836          ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
 837         else
 838          sslipi=1.0d0
 839          ssgradlipi=0.0
 840         endif
 841        else
 842          sslipi=0.0d0
 843          ssgradlipi=0.0
 844        endif
 845
 846         dxi=dc_norm(1,nres+i)
 847         dyi=dc_norm(2,nres+i)
 848         dzi=dc_norm(3,nres+i)
 849         dsci_inv=dsc_inv(itypi)
 850 C
 851 C Calculate SC interaction energy.
 852 C
 853 c        do iint=1,nint_gr(i)
 854 c          do j=istart(i,iint),iend(i,iint)
 855          do j=i+1,nct
 856           IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
 857             ind=ind+1
 858             itypj=iabs(itype(j))
 859             if (itypj.eq.ntyp1) cycle
 860             dscj_inv=dsc_inv(itypj)
 861             sig0ij=sigma(itypi,itypj)
 862             chi1=chi(itypi,itypj)
 863             chi2=chi(itypj,itypi)
 864             chi12=chi1*chi2
 865             chip1=chip(itypi)
 866             chip2=chip(itypj)
 867             chip12=chip1*chip2
 868             alf1=alp(itypi)
 869             alf2=alp(itypj)
 870             alf12=0.5D0*(alf1+alf2)
 871 C For diagnostics only!!!
 872 c           chi1=0.0D0
 873 c           chi2=0.0D0
 874 c           chi12=0.0D0
 875 c           chip1=0.0D0
 876 c           chip2=0.0D0
 877 c           chip12=0.0D0
 878 c           alf1=0.0D0
 879 c           alf2=0.0D0
 880 c           alf12=0.0D0
 881 C            xj=c(1,nres+j)-xi
 882 C            yj=c(2,nres+j)-yi
 883 C            zj=c(3,nres+j)-zi
 884             xj=c(1,nres+j)
 885             yj=c(2,nres+j)
 886             zj=c(3,nres+j)
 887           xj=mod(xj,boxxsize)
 888           if (xj.lt.0) xj=xj+boxxsize
 889           yj=mod(yj,boxysize)
 890           if (yj.lt.0) yj=yj+boxysize
 891           zj=mod(zj,boxzsize)
 892           if (zj.lt.0) zj=zj+boxzsize
 893        if ((zj.gt.bordlipbot)
 894      &.and.(zj.lt.bordliptop)) then
 895 C the energy transfer exist
 896         if (zj.lt.buflipbot) then
 897 C what fraction I am in
 898          fracinbuf=1.0d0-
 899      &        ((positi-bordlipbot)/lipbufthick)
 900 C lipbufthick is thickenes of lipid buffore
 901          sslipj=sscalelip(fracinbuf)
 902          ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
 903         elseif (zi.gt.bufliptop) then
 904          fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
 905          sslipj=sscalelip(fracinbuf)
 906          ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
 907         else
 908          sslipj=1.0d0
 909          ssgradlipj=0.0
 910         endif
 911        else
 912          sslipj=0.0d0
 913          ssgradlipj=0.0
 914        endif
 915       aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
 916      &  +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
 917       bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
 918      &  +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
 919
 920       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
 921       xj_safe=xj
 922       yj_safe=yj
 923       zj_safe=zj
 924       subchap=0
 925       do xshift=-1,1
 926       do yshift=-1,1
 927       do zshift=-1,1
 928           xj=xj_safe+xshift*boxxsize
 929           yj=yj_safe+yshift*boxysize
 930           zj=zj_safe+zshift*boxzsize
 931           dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
 932           if(dist_temp.lt.dist_init) then
 933             dist_init=dist_temp
 934             xj_temp=xj
 935             yj_temp=yj
 936             zj_temp=zj
 937             subchap=1
 938           endif
 939        enddo
 940        enddo
 941        enddo
 942        if (subchap.eq.1) then
 943           xj=xj_temp-xi
 944           yj=yj_temp-yi
 945           zj=zj_temp-zi
 946        else
 947           xj=xj_safe-xi
 948           yj=yj_safe-yi
 949           zj=zj_safe-zi
 950        endif
 951
 952             dxj=dc_norm(1,nres+j)
 953             dyj=dc_norm(2,nres+j)
 954             dzj=dc_norm(3,nres+j)
 955             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
 956             rij=dsqrt(rrij)
 957             sss=sscale((1.0d0/rij)/sigma(itypi,itypj))
 958             sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj))
 959
 960 C Calculate angle-dependent terms of energy and contributions to their
 961 C derivatives.
 962             call sc_angular
 963             sigsq=1.0D0/sigsq
 964             sig=sig0ij*dsqrt(sigsq)
 965             rij_shift=1.0D0/rij-sig+sig0ij
 966 C I hate to put IF's in the loops, but here don't have another choice!!!!
 967             if (rij_shift.le.0.0D0) then
 968               evdw=1.0D20
 969 cd              write (iout,'(2(a3,i3,2x),17(0pf7.3))')
 970 cd     &        restyp(itypi),i,restyp(itypj),j,
 971 cd     &        rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
 972               return
 973             endif
 974             sigder=-sig*sigsq
 975 c---------------------------------------------------------------
 976             rij_shift=1.0D0/rij_shift
 977             fac=rij_shift**expon
 978             e1=fac*fac*aa
 979             e2=fac*bb
 980             evdwij=eps1*eps2rt*eps3rt*(e1+e2)
 981             eps2der=evdwij*eps3rt
 982             eps3der=evdwij*eps2rt
 983             evdwij=evdwij*eps2rt*eps3rt
 984             evdw=evdw+evdwij
 985             if (lprn) then
 986             sigm=dabs(aa/bb)**(1.0D0/6.0D0)
 987             epsi=bb**2/aa
 988 cd            write (iout,'(2(a3,i3,2x),17(0pf7.3))')
 989 cd     &        restyp(itypi),i,restyp(itypj),j,
 990 cd     &        epsi,sigm,chi1,chi2,chip1,chip2,
 991 cd     &        eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
 992 cd     &        om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
 993 cd     &        evdwij
 994             endif
 995
 996             if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
 997      &                        'evdw',i,j,evdwij
 998
 999 C Calculate gradient components.
1000             e1=e1*eps1*eps2rt**2*eps3rt**2
1001             fac=-expon*(e1+evdwij)*rij_shift
1002             sigder=fac*sigder
1003             fac=rij*fac
1004             fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij
1005 C Calculate the radial part of the gradient
1006             gg(1)=xj*fac
1007             gg(2)=yj*fac
1008             gg(3)=zj*fac
1009             gg_lipi(3)=ssgradlipi*evdwij
1010             gg_lipj(3)=ssgradlipj*evdwij
1011 C Calculate angular part of the gradient.
1012             call sc_grad
1013           ENDIF
1014           enddo      ! j
1015 c        enddo        ! iint
1016       enddo          ! i
1017       end
1018 C-----------------------------------------------------------------------------