source/unres/src_CSA/minim_jlee.F

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