source/unres/src_CSA_DiL/minim_jlee.F

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