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