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