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