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