source/unres/src-HCD-5D/minim_jlee.F

   1       subroutine minim_jlee
   2 #ifdef LBFGS
   3       use minima
   4       use inform
   5       use output
   6       use iounit
   7       use scales
   8 #endif
   9 c  controls minimization and sorting routines
  10       implicit real*8 (a-h,o-z)
  11       include 'DIMENSIONS'
  12 #ifndef LBFGS
  13       integer liv,lv
  14       parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
  15 #endif
  16       include 'COMMON.VAR'
  17       include 'COMMON.IOUNITS'
  18       include 'COMMON.MINIM'
  19       include 'COMMON.CONTROL'
  20 #ifdef LBFGS
  21       common /gacia/ nfun
  22       double precision grdmin
  23       external funcgrad
  24       external optsave
  25 #else
  26       external func,gradient,fdum
  27       dimension iv(liv)
  28       double precision v(1:lv+1)
  29       common /przechowalnia/ v
  30 #endif
  31       real ran1,ran2,ran3
  32 #ifdef MPI
  33       include 'mpif.h'
  34       include 'COMMON.SETUP'
  35       dimension muster(mpi_status_size)
  36 #endif
  37       include 'COMMON.GEO'
  38       include 'COMMON.FFIELD'
  39       include 'COMMON.SBRIDGE'
  40       include 'COMMON.DISTFIT'
  41       include 'COMMON.CHAIN'
  42       dimension var(maxvar),erg(mxch*(mxch+1)/2+1)
  43       dimension var2(maxvar)
  44       integer iffr(maxres),ihpbt(maxdim),jhpbt(maxdim)
  45       double precision d(maxvar),garbage(maxvar),g(maxvar)
  46       double precision energia(0:n_ene),time0s,time1s
  47       dimension indx(9),info(12)
  48       dimension idum(1),rdum(1)
  49       dimension icont(2,maxcont)
  50       logical check_var,fail
  51       integer iloop(2)
  52       data rad /1.745329252d-2/
  53 c  receive # of start
  54 !      print *,'Processor',me,' calling MINIM_JLEE maxfun',maxfun,
  55 !     &   ' maxmin',maxmin,' tolf',tolf,' rtolf',rtolf
  56 #ifdef LBFGS
  57       maxiter=maxmin
  58       coordtype='RIGIDBODY'
  59       grdmin=tolf
  60       jout=iout
  61       jprint=print_min_stat
  62       iwrite=0
  63       if (.not. allocated(scale))  allocate (scale(nvar))
  64 c
  65 c     set scaling parameter for function and derivative values;
  66 c     use square root of median eigenvalue of typical Hessian
  67 c
  68       set_scale = .true.
  69 c      nvar = 0
  70       do i = 1, nvar
  71 c         if (use(i)) then
  72 c            do j = 1, 3
  73 c               nvar = nvar + 1
  74                scale(i) = 12.0d0
  75 c            end do
  76 c         end if
  77       end do
  78 #endif
  79       nhpb0=nhpb
  80    10 continue
  81       time0s=MPI_WTIME()
  82 c      print *, 'MINIM_JLEE: ',me,' is waiting'
  83       call mpi_recv(info,12,mpi_integer,king,idint,CG_COMM,
  84      *              muster,ierr)
  85       time1s=MPI_WTIME()
  86       write (iout,'(a12,f10.4,a4)')'Waiting for ',time1s-time0s,' sec'
  87       call flush(iout)
  88        n=info(1)
  89 c      print *, 'MINIM_JLEE: ',me,' received: ',n
  90
  91 crc      if (ierr.ne.0) go to 100
  92 c  if # = 0, return
  93       if (n.eq.0) then
  94         write (iout,*) 'Finishing minim_jlee - signal',n,' from master'
  95         call flush(iout)
  96         return
  97       endif
  98
  99       nfun=0
 100       IF (n.lt.0) THEN
 101        call mpi_recv(var,nvar,mpi_double_precision,
 102      *              king,idreal,CG_COMM,muster,ierr)
 103        call mpi_recv(iffr,nres,mpi_integer,
 104      *              king,idint,CG_COMM,muster,ierr)
 105        call mpi_recv(var2,nvar,mpi_double_precision,
 106      *              king,idreal,CG_COMM,muster,ierr)
 107       ELSE
 108 c  receive initial values of variables
 109        call mpi_recv(var,nvar,mpi_double_precision,
 110      *              king,idreal,CG_COMM,muster,ierr)
 111 crc       if (ierr.ne.0) go to 100
 112       ENDIF
 113
 114       if(vdisulf.and.info(2).ne.-1) then
 115         if(info(4).ne.0)then
 116            call mpi_recv(ihpbt,info(4),mpi_integer,
 117      *              king,idint,CG_COMM,muster,ierr)
 118            call mpi_recv(jhpbt,info(4),mpi_integer,
 119      *              king,idint,CG_COMM,muster,ierr)
 120         endif
 121       endif
 122
 123       IF (n.lt.0) THEN
 124         n=-n
 125         nhpb=nhpb0
 126         link_start=1
 127         link_end=nhpb
 128         call init_int_table
 129         call contact_cp(var,var2,iffr,nfun,n)
 130       ENDIF
 131
 132       if(vdisulf.and.info(2).ne.-1) then
 133         nss=0
 134         if(info(4).ne.0)then
 135 cd          write(iout,*) 'SS=',info(4),'N=',info(1),'IT=',info(2)
 136           call var_to_geom(nvar,var)
 137           call chainbuild
 138           do i=1,info(4)
 139            if (dist(ihpbt(i),jhpbt(i)).lt.7.0) then
 140             nss=nss+1
 141             ihpb(nss)=ihpbt(i)
 142             jhpb(nss)=jhpbt(i)
 143 cd            write(iout,*) 'SS  mv=',info(3),
 144 cd     &         ihpb(nss)-nres,jhpb(nss)-nres,
 145 cd     &         dist(ihpb(nss),jhpb(nss))
 146             dhpb(nss)=dbr
 147             forcon(nss)=fbr
 148            else
 149 cd            write(iout,*) 'rm SS  mv=',info(3),
 150 cd     &         ihpbt(i)-nres,jhpbt(i)-nres,dist(ihpbt(i),jhpbt(i))
 151            endif
 152           enddo
 153         endif
 154         nhpb=nss
 155         link_start=1
 156         link_end=nhpb
 157         call init_int_table
 158       endif
 159
 160       if (info(3).eq.14) then
 161        write(iout,*) 'calling local_move',info(7),info(8)
 162        call local_move_init(.false.)
 163        call var_to_geom(nvar,var)
 164        call local_move(info(7),info(8),20d0,50d0)
 165        call geom_to_var(nvar,var)
 166       endif
 167
 168
 169       if (info(3).eq.16) then
 170        write(iout,*) 'calling beta_slide',info(7),info(8),
 171      &            info(10), info(11), info(12)
 172        call var_to_geom(nvar,var)
 173        call beta_slide(info(7),info(8),info(10),info(11),info(12)
 174      &                                                     ,nfun,n)
 175        call geom_to_var(nvar,var)
 176       endif
 177
 178
 179       if (info(3).eq.17) then
 180        write(iout,*) 'calling beta_zip',info(7),info(8)
 181        call var_to_geom(nvar,var)
 182        call beta_zip(info(7),info(8),nfun,n)
 183        call geom_to_var(nvar,var)
 184       endif
 185
 186
 187 crc overlap test
 188
 189       if (overlapsc) then
 190
 191           call var_to_geom(nvar,var)
 192           call chainbuild
 193           call etotal(energia(0))
 194           nfun=nfun+1
 195           if (energia(1).eq.1.0d20) then
 196             info(3)=-info(3)
 197             write (iout,'(a,1pe14.5)')'#OVERLAP evdw=1d20',energia(1)
 198             call overlap_sc(fail)
 199             if(.not.fail) then
 200              call geom_to_var(nvar,var)
 201              call etotal(energia(0))
 202              nfun=nfun+1
 203              write (iout,'(a,1pe14.5)')'#OVERLAP evdw after',energia(1)
 204             else
 205 #ifdef LBFGS
 206              etot=1.0d20
 207              nfun=-1
 208 #else
 209              v(10)=1.0d20
 210              iv(1)=-1
 211 #endif
 212              goto 201
 213             endif
 214           endif
 215       endif
 216
 217       if (searchsc) then
 218           call var_to_geom(nvar,var)
 219           call sc_move(2,nres-1,1,10d0,nft_sc,etot)
 220           call geom_to_var(nvar,var)
 221 cd          write(iout,*) 'sc_move',nft_sc,etot
 222       endif
 223
 224       if (check_var(var,info)) then
 225 #ifdef LBFGS
 226            etot=1.0d21
 227 #else
 228            v(10)=1.0d21
 229            iv(1)=6
 230 #endif
 231            goto 201
 232       endif
 233
 234
 235 crc
 236
 237 !      write (iout,*) 'MINIM_JLEE: Processor',me,' nvar',nvar
 238 !      write (iout,'(8f10.4)') (var(i),i=1,nvar)
 239 !      write (*,*) 'MINIM_JLEE: Processor',me,' received nvar',nvar
 240 !      write (*,'(8f10.4)') (var(i),i=1,nvar)
 241
 242       do i=1,nvar
 243         garbage(i)=var(i)
 244       enddo
 245 #ifdef LBFGS
 246       eee=funcgrad(var,g)
 247       nfun=nfun+1
 248       if(eee.ge.1.0d20) then
 249 c       print *,'MINIM_JLEE: ',me,' CHUJ NASTAPIL'
 250 c       print *,' energy before SUMSL =',eee
 251 c       print *,' aborting local minimization'
 252        go to 201
 253       endif
 254
 255       call lbfgs (nvar,var,etot,grdmin,funcgrad,optsave)
 256       deallocate(scale)
 257 #else
 258       call deflt(2,iv,liv,lv,v)
 259 * 12 means fresh start, dont call deflt
 260       iv(1)=12
 261 * max num of fun calls
 262       if (maxfun.eq.0) maxfun=500
 263       iv(17)=maxfun
 264 * max num of iterations
 265       if (maxmin.eq.0) maxmin=1000
 266       iv(18)=maxmin
 267 * controls output
 268       iv(19)=2
 269 * selects output unit
 270 cd      iv(21)=iout
 271       iv(21)=0
 272 * 1 means to print out result
 273       iv(22)=0
 274 cd        iv(22)=1
 275 * 1 means to print out summary stats
 276       iv(23)=0
 277 * 1 means to print initial x and d
 278       iv(24)=0
 279
 280 c      if(me.eq.3.and.n.eq.255) then
 281 c       print *,' CHUJ: stoi'
 282 c       iv(21)=6
 283 c       iv(22)=1
 284 c       iv(23)=1
 285 c       iv(24)=1
 286 c      endif
 287
 288 * min val for v(radfac) default is 0.1
 289       v(24)=0.1D0
 290 * max val for v(radfac) default is 4.0
 291       v(25)=2.0D0
 292 c     v(25)=4.0D0
 293 * check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
 294 * the sumsl default is 0.1
 295       v(26)=0.1D0
 296 * false conv if (act fnctn decrease) .lt. v(34)
 297 * the sumsl default is 100*machep
 298       v(34)=v(34)/100.0D0
 299 * absolute convergence
 300       if (tolf.eq.0.0D0) tolf=1.0D-4
 301       v(31)=tolf
 302 * relative convergence
 303       if (rtolf.eq.0.0D0) rtolf=1.0D-4
 304       v(32)=rtolf
 305 * controls initial step size
 306        v(35)=1.0D-1
 307 * large vals of d correspond to small components of step
 308       do i=1,nphi
 309         d(i)=1.0D-1
 310       enddo
 311       do i=nphi+1,nvar
 312         d(i)=1.0D-1
 313       enddo
 314 c  minimize energy
 315 !      write (iout,*) 'Processor',me,' nvar',nvar
 316 !      write (iout,*) 'Variables BEFORE minimization:'
 317 !      write (iout,'(8f10.4)') (rad2deg*var(i),i=1,nvar)
 318
 319 c      print *, 'MINIM_JLEE: ',me,' before SUMSL '
 320
 321       call func(nvar,var,nf,eee,idum,rdum,fdum)
 322       nfun=nfun+1
 323       if(eee.ge.1.0d20) then
 324 c       print *,'MINIM_JLEE: ',me,' CHUJ NASTAPIL'
 325 c       print *,' energy before SUMSL =',eee
 326 c       print *,' aborting local minimization'
 327 #ifdef LBFGS
 328        etot=eee
 329 #else
 330        iv(1)=-1
 331        v(10)=eee
 332 #endif
 333        go to 201
 334       endif
 335
 336 ct      time0s=MPI_WTIME()
 337       call sumsl(nvar,d,var,func,gradient,iv,liv,lv,v,idum,rdum,fdum)
 338 ct      write(iout,*) 'sumsl time=',MPI_WTIME()-time0s,iv(7),v(10)
 339 c      print *, 'MINIM_JLEE: ',me,' after SUMSL '
 340
 341 c  find which conformation was returned from sumsl
 342         nfun=nfun+iv(7)
 343 #endif
 344 !      print *,'Processor',me,' iv(17)',iv(17),' iv(18)',iv(18),' nf',nf,
 345 !     & ' retcode',iv(1),' energy',v(10),' tolf',v(31),' rtolf',v(32)
 346 c        if (iv(1).ne.4 .or. nf.le.1) then
 347 c          write (*,*) 'Processor',me,' something bad in SUMSL',iv(1),nf
 348 c         write (*,*) 'Initial Variables'
 349 c         write (*,'(8f10.4)') (rad2deg*garbage(i),i=1,nvar)
 350 c         write (*,*) 'Variables'
 351 c         write (*,'(8f10.4)') (rad2deg*var(i),i=1,nvar)
 352 c         write (*,*) 'Vector d'
 353 c         write (*,'(8f10.4)') (d(i),i=1,nvar)
 354 c         write (iout,*) 'Processor',me,' something bad in SUMSL',
 355 c    &       iv(1),nf
 356 c         write (iout,*) 'Initial Variables'
 357 c         write (iout,'(8f10.4)') (rad2deg*garbage(i),i=1,nvar)
 358 c         write (iout,*) 'Variables'
 359 c         write (iout,'(8f10.4)') (rad2deg*var(i),i=1,nvar)
 360 c         write (iout,*) 'Vector d'
 361 c         write (iout,'(8f10.4)') (d(i),i=1,nvar)
 362 c        endif
 363 c        if (nf.lt.iv(6)-1) then
 364 c  recalculate intra- and interchain energies
 365 c         call func(nvar,var,nf,v(10),iv,v,fdum)
 366 c        else if (nf.eq.iv(6)-1) then
 367 c  regenerate conformation
 368 c         call var_to_geom(nvar,var)
 369 c         call chainbuild
 370 c        endif
 371 c  change origin and axes to standard ECEPP format
 372 c        call var_to_geom(nvar,var)
 373 !      write (iout,*) 'MINIM_JLEE after minim: Processor',me,' nvar',nvar
 374 !      write (iout,'(8f10.4)') (var(i),i=1,nvar)
 375 !      write (iout,*) 'Energy:',v(10)
 376 c  send back output
 377 c       print *, 'MINIM_JLEE: ',me,' minimized: ',n
 378   201  continue
 379         indx(1)=n
 380 c return code: 6-gradient 9-number of ftn evaluation, etc
 381 #ifdef LBFGS
 382         indx(2)=nfun
 383 #else
 384         indx(2)=iv(1)
 385 #endif
 386 c total # of ftn evaluations (for iwf=0, it includes all minimizations).
 387         indx(3)=nfun
 388         indx(4)=info(2)
 389         indx(5)=info(3)
 390         indx(6)=nss
 391         indx(7)=info(5)
 392         indx(8)=info(6)
 393         indx(9)=info(9)
 394         call mpi_send(indx,9,mpi_integer,king,idint,CG_COMM,
 395      *                 ierr)
 396 c  send back energies
 397 c al & cc
 398 c calculate contact order
 399 #ifdef LBFGS
 400 #ifdef CO_BIAS
 401         call contact(.false.,ncont,icont,co)
 402         erg(1)=etot-1.0d2*co
 403 #else
 404         erg(1)=etot
 405 #endif
 406 #else
 407 #ifdef CO_BIAS
 408         call contact(.false.,ncont,icont,co)
 409         erg(1)=v(10)-1.0d2*co
 410 #else
 411         erg(1)=v(10)
 412 #endif
 413 #endif
 414         j=1
 415         call mpi_send(erg,j,mpi_double_precision,king,idreal,
 416      *                 CG_COMM,ierr)
 417 #ifdef CO_BIAS
 418         call mpi_send(co,j,mpi_double_precision,king,idreal,
 419      *                 CG_COMM,ierr)
 420 #endif
 421 c  send back values of variables
 422         call mpi_send(var,nvar,mpi_double_precision,
 423      *               king,idreal,CG_COMM,ierr)
 424 !        print * , 'MINIM_JLEE: Processor',me,' send erg and var '
 425
 426         if(vdisulf.and.info(2).ne.-1.and.nss.ne.0) then
 427 cd          call intout
 428 cd          call chainbuild
 429 cd          call etotal(energia(0))
 430 cd          etot=energia(0)
 431 cd          call enerprint(energia(0))
 432          call mpi_send(ihpb,nss,mpi_integer,
 433      *               king,idint,CG_COMM,ierr)
 434          call mpi_send(jhpb,nss,mpi_integer,
 435      *               king,idint,CG_COMM,ierr)
 436         endif
 437
 438         go to 10
 439   100 print *, ' error in receiving message from emperor', me
 440       call mpi_abort(mpi_comm_world,ierror,ierrcode)
 441       return
 442   200 print *, ' error in sending message to emperor'
 443       call mpi_abort(mpi_comm_world,ierror,ierrcode)
 444       return
 445   300 print *, ' error in communicating with emperor'
 446       call mpi_abort(mpi_comm_world,ierror,ierrcode)
 447       return
 448   956 format (' initial energy could not be calculated',41x)
 449   957 format (80x)
 450   965 format (' convergence code ',i2,' # of function calls ',
 451      *  i4,' # of gradient calls ',i4,10x)
 452   975 format (' energy ',1p,e12.4,' scaled gradient ',e11.3,32x)
 453       end
 454
 455       logical function check_var(var,info)
 456       implicit real*8 (a-h,o-z)
 457       include 'DIMENSIONS'
 458       include 'COMMON.VAR'
 459       include 'COMMON.IOUNITS'
 460       include 'COMMON.GEO'
 461       include 'COMMON.SETUP'
 462       dimension var(maxvar)
 463       dimension info(3)
 464 C AL -------
 465        check_var=.false.
 466        do i=nphi+ntheta+1,nphi+ntheta+nside
 467 ! Check the side chain "valence" angles alpha
 468          if (var(i).lt.1.0d-7) then
 469            write (iout,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
 470            write (iout,*) 'Processor',me,'received bad variables!!!!'
 471            write (iout,*) 'Variables'
 472            write (iout,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
 473            write (iout,*) 'Continuing calculations at this point',
 474      & ' could destroy the results obtained so far... ABORTING!!!!!!'
 475            write (iout,'(a19,i5,f10.4,a4,2i4,a3,i3)')
 476      &     'valence angle alpha',i-nphi-ntheta,var(i),
 477      &     'n it',info(1),info(2),'mv ',info(3)
 478            write (*,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
 479            write (*,*) 'Processor',me,'received bad variables!!!!'
 480            write (*,*) 'Variables'
 481            write (*,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
 482            write (*,*) 'Continuing calculations at this point',
 483      & ' could destroy the results obtained so far... ABORTING!!!!!!'
 484            write (*,'(a19,i5,f10.4,a4,2i4,a3,i3)')
 485      &     'valence angle alpha',i-nphi-ntheta,var(i),
 486      &     'n it',info(1),info(2),'mv ',info(3)
 487            check_var=.true.
 488            return
 489          endif
 490        enddo
 491 ! Check the backbone "valence" angles theta
 492        do i=nphi+1,nphi+ntheta
 493          if (var(i).lt.1.0d-7) then
 494            write (iout,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
 495            write (iout,*) 'Processor',me,'received bad variables!!!!'
 496            write (iout,*) 'Variables'
 497            write (iout,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
 498            write (iout,*) 'Continuing calculations at this point',
 499      & ' could destroy the results obtained so far... ABORTING!!!!!!'
 500            write (iout,'(a19,i5,f10.4,a4,2i4,a3,i3)')
 501      &     'valence angle theta',i-nphi,var(i),
 502      &     'n it',info(1),info(2),'mv ',info(3)
 503            write (*,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
 504            write (*,*) 'Processor',me,'received bad variables!!!!'
 505            write (*,*) 'Variables'
 506            write (*,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
 507            write (*,*) 'Continuing calculations at this point',
 508      & ' could destroy the results obtained so far... ABORTING!!!!!!'
 509            write (*,'(a19,i5,f10.4,a4,2i4,a3,i3)')
 510      &     'valence angle theta',i-nphi,var(i),
 511      &     'n it',info(1),info(2),'mv ',info(3)
 512            check_var=.true.
 513            return
 514          endif
 515        enddo
 516        return
 517        end