source/unres/src-HCD-5D/minimize_p.F

   1       subroutine minimize(etot,x,iretcode,nfun)
   2       implicit none
   3       include 'DIMENSIONS'
   4       integer liv,lv
   5       parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
   6 *********************************************************************
   7 * OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
   8 * the calling subprogram.                                           *
   9 * when d(i)=1.0, then v(35) is the length of the initial step,      *
  10 * calculated in the usual pythagorean way.                          *
  11 * absolute convergence occurs when the function is within v(31) of  *
  12 * zero. unless you know the minimum value in advance, abs convg     *
  13 * is probably not useful.                                           *
  14 * relative convergence is when the model predicts that the function *
  15 * will decrease by less than v(32)*abs(fun).                        *
  16 *********************************************************************
  17       include 'COMMON.IOUNITS'
  18       include 'COMMON.VAR'
  19       include 'COMMON.GEO'
  20       include 'COMMON.MINIM'
  21       integer icall
  22       common /srutu/ icall
  23       integer iv(liv)
  24       double precision minval,x(maxvar),d(maxvar),v(1:lv),xx(maxvar)
  25       double precision energia(0:n_ene)
  26       integer idum
  27       double precision rdum
  28       double precision fdum
  29       external func,gradient,fdum
  30       external func_restr,grad_restr
  31       logical not_done,change,reduce
  32       integer i,nvar_restr,nfun,iretcode
  33       double precision etot
  34 c      common /przechowalnia/ v
  35
  36       icall = 1
  37
  38       NOT_DONE=.TRUE.
  39
  40 c     DO WHILE (NOT_DONE)
  41
  42       call deflt(2,iv,liv,lv,v)
  43 * 12 means fresh start, dont call deflt
  44       iv(1)=12
  45 * max num of fun calls
  46       if (maxfun.eq.0) maxfun=500
  47       iv(17)=maxfun
  48 * max num of iterations
  49       if (maxmin.eq.0) maxmin=1000
  50       iv(18)=maxmin
  51 * controls output
  52       iv(19)=2
  53 * selects output unit
  54       iv(21)=0
  55       if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout
  56 * 1 means to print out result
  57       iv(22)=print_min_res
  58 * 1 means to print out summary stats
  59       iv(23)=print_min_stat
  60 * 1 means to print initial x and d
  61       iv(24)=print_min_ini
  62 * min val for v(radfac) default is 0.1
  63       v(24)=0.1D0
  64 * max val for v(radfac) default is 4.0
  65       v(25)=2.0D0
  66 c     v(25)=4.0D0
  67 * check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
  68 * the sumsl default is 0.1
  69       v(26)=0.1D0
  70 * false conv if (act fnctn decrease) .lt. v(34)
  71 * the sumsl default is 100*machep
  72       v(34)=v(34)/100.0D0
  73 * absolute convergence
  74       if (tolf.eq.0.0D0) tolf=1.0D-4
  75       v(31)=tolf
  76 * relative convergence
  77       if (rtolf.eq.0.0D0) rtolf=1.0D-4
  78       v(32)=rtolf
  79 * controls initial step size
  80       v(35)=1.0D-1
  81 * large vals of d correspond to small components of step
  82       do i=1,nphi
  83         d(i)=1.0D-1
  84       enddo
  85       do i=nphi+1,nvar
  86         d(i)=1.0D-1
  87       enddo
  88 cd    print *,'Calling SUMSL'
  89 c     call var_to_geom(nvar,x)
  90 c     call chainbuild
  91 c     call etotal(energia(0))
  92 c     etot = energia(0)
  93       IF (mask_r) THEN
  94        call x2xx(x,xx,nvar_restr)
  95        call sumsl(nvar_restr,d,xx,func_restr,grad_restr,
  96      &                    iv,liv,lv,v,idum,rdum,fdum)
  97        call xx2x(x,xx)
  98       ELSE
  99        call sumsl(nvar,d,x,func,gradient,iv,liv,lv,v,idum,rdum,fdum)
 100       ENDIF
 101       etot=v(10)
 102       iretcode=iv(1)
 103 cd    print *,'Exit SUMSL; return code:',iretcode,' energy:',etot
 104 cd    write (iout,'(/a,i4/)') 'SUMSL return code:',iv(1)
 105 c     call intout
 106 c     change=reduce(x)
 107       call var_to_geom(nvar,x)
 108 c     if (change) then
 109 c       write (iout,'(a)') 'Reduction worked, minimizing again...'
 110 c     else
 111 c       not_done=.false.
 112 c     endif
 113       call chainbuild
 114 c     call etotal(energia(0))
 115 c     etot=energia(0)
 116 c     call enerprint(energia(0))
 117       nfun=iv(6)
 118
 119 c     write (*,*) 'Processor',MyID,' leaves MINIMIZE.'
 120
 121 c     ENDDO ! NOT_DONE
 122
 123       return
 124       end
 125 #ifdef MPI
 126 c----------------------------------------------------------------------------
 127       subroutine ergastulum
 128       implicit none
 129       include 'DIMENSIONS'
 130 #ifdef MPI
 131       include "mpif.h"
 132       double precision time00
 133       integer ierr,ierror
 134 #endif
 135       include 'COMMON.SETUP'
 136       include 'COMMON.DERIV'
 137       include 'COMMON.VAR'
 138       include 'COMMON.IOUNITS'
 139       include 'COMMON.FFIELD'
 140       include 'COMMON.INTERACT'
 141       include 'COMMON.MD'
 142 #ifdef FIVEDIAG
 143        include 'COMMON.LAGRANGE.5diag'
 144 #else
 145        include 'COMMON.LAGRANGE'
 146 #endif
 147       include 'COMMON.TIME1'
 148       double precision z(maxres6),d_a_tmp(maxres6)
 149       double precision edum(0:n_ene),time_order(0:10)
 150       double precision Gcopy(maxres2,maxres2)
 151       common /przechowalnia/ Gcopy
 152       integer icall /0/
 153       integer i,j,iorder
 154 C Workers wait for variables and NF, and NFL from the boss
 155       iorder=0
 156       do while (iorder.ge.0)
 157 c      write (*,*) 'Processor',fg_rank,' CG group',kolor,
 158 c     & ' receives order from Master'
 159 #ifdef MPI
 160         time00=MPI_Wtime()
 161         call MPI_Bcast(iorder,1,MPI_INTEGER,king,FG_COMM,IERR)
 162         time_Bcast=time_Bcast+MPI_Wtime()-time00
 163         if (icall.gt.4 .and. iorder.ge.0)
 164      &   time_order(iorder)=time_order(iorder)+MPI_Wtime()-time00
 165 #endif
 166        icall=icall+1
 167 c      write (*,*)
 168 c     & 'Processor',fg_rank,' completed receive MPI_BCAST order',iorder
 169         if (iorder.eq.0) then
 170           call zerograd
 171           call etotal(edum)
 172 c          write (2,*) "After etotal"
 173 c          write (2,*) "dimen",dimen," dimen3",dimen3
 174 c          call flush(2)
 175         else if (iorder.eq.2) then
 176           call zerograd
 177           call etotal_short(edum)
 178 c          write (2,*) "After etotal_short"
 179 c          write (2,*) "dimen",dimen," dimen3",dimen3
 180 c          call flush(2)
 181         else if (iorder.eq.3) then
 182           call zerograd
 183           call etotal_long(edum)
 184 c          write (2,*) "After etotal_long"
 185 c          write (2,*) "dimen",dimen," dimen3",dimen3
 186 c          call flush(2)
 187         else if (iorder.eq.1) then
 188           call sum_gradient
 189 c          write (2,*) "After sum_gradient"
 190 c          write (2,*) "dimen",dimen," dimen3",dimen3
 191 c          call flush(2
 192 #ifndef FIVEDIAG
 193         else if (iorder.eq.4) then
 194           call ginv_mult(z,d_a_tmp)
 195         else if (iorder.eq.5) then
 196 c Setup MD things for a slave
 197           dimen=(nct-nnt+1)+nside
 198           dimen1=(nct-nnt)+(nct-nnt+1)
 199           dimen3=dimen*3
 200 c          write (2,*) "dimen",dimen," dimen3",dimen3
 201 c          call flush(2)
 202           call int_bounds(dimen,igmult_start,igmult_end)
 203           igmult_start=igmult_start-1
 204           call MPI_Allgather(3*igmult_start,1,MPI_INTEGER,
 205      &     ng_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
 206            my_ng_count=igmult_end-igmult_start
 207           call MPI_Allgather(3*my_ng_count,1,MPI_INTEGER,ng_counts(0),1,
 208      &     MPI_INTEGER,FG_COMM,IERROR)
 209 c          write (2,*) "ng_start",(ng_start(i),i=0,nfgtasks-1)
 210 c          write (2,*) "ng_counts",(ng_counts(i),i=0,nfgtasks-1)
 211           myginv_ng_count=maxres2*my_ng_count
 212 c          write (2,*) "igmult_start",igmult_start," igmult_end",
 213 c     &     igmult_end," my_ng_count",my_ng_count
 214 c          call flush(2)
 215           call MPI_Allgather(maxres2*igmult_start,1,MPI_INTEGER,
 216      &     nginv_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
 217           call MPI_Allgather(myginv_ng_count,1,MPI_INTEGER,
 218      &     nginv_counts(0),1,MPI_INTEGER,FG_COMM,IERROR)
 219 c          write (2,*) "nginv_start",(nginv_start(i),i=0,nfgtasks-1)
 220 c          write (2,*) "nginv_counts",(nginv_counts(i),i=0,nfgtasks-1)
 221 c          call flush(2)
 222 c          call MPI_Barrier(FG_COMM,IERROR)
 223           time00=MPI_Wtime()
 224           call MPI_Scatterv(ginv(1,1),nginv_counts(0),
 225      &     nginv_start(0),MPI_DOUBLE_PRECISION,gcopy(1,1),
 226      &     myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 227 #ifdef TIMING
 228           time_scatter_ginv=time_scatter_ginv+MPI_Wtime()-time00
 229 #endif
 230           do i=1,dimen
 231             do j=1,2*my_ng_count
 232               ginv(j,i)=gcopy(i,j)
 233             enddo
 234           enddo
 235 c          write (2,*) "dimen",dimen," dimen3",dimen3
 236 c          write (2,*) "End MD setup"
 237 c          call flush(2)
 238 c           write (iout,*) "My chunk of ginv_block"
 239 c           call MATOUT2(my_ng_count,dimen3,maxres2,maxers2,ginv_block)
 240 #endif
 241         else if (iorder.eq.6) then
 242           call int_from_cart1(.false.)
 243         else if (iorder.eq.7) then
 244           call chainbuild_cart
 245         else if (iorder.eq.8) then
 246           call intcartderiv
 247 #ifndef FIVEDIAG
 248         else if (iorder.eq.9) then
 249           call fricmat_mult(z,d_a_tmp)
 250 #endif
 251         else if (iorder.eq.10) then
 252           call setup_fricmat
 253         endif
 254       enddo
 255       write (*,*) 'Processor',fg_rank,' CG group',kolor,
 256      &  ' absolute rank',myrank,' leves ERGASTULUM.'
 257       write(*,*)'Processor',fg_rank,' wait times for respective orders',
 258      & (' order[',i,']',time_order(i),i=0,10)
 259       return
 260       end
 261 #endif
 262 ************************************************************************
 263       subroutine func(n,x,nf,f,uiparm,urparm,ufparm)
 264       implicit real*8 (a-h,o-z)
 265       include 'DIMENSIONS'
 266       include 'COMMON.DERIV'
 267       include 'COMMON.IOUNITS'
 268       include 'COMMON.GEO'
 269       common /chuju/ jjj
 270       double precision energia(0:n_ene)
 271       integer jjj
 272       double precision ufparm
 273       external ufparm
 274       integer uiparm(1)
 275       real*8 urparm(1)
 276       dimension x(maxvar)
 277 c     if (jjj.gt.0) then
 278 c      write (iout,*) "in func x"
 279 c      write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
 280 c     endif
 281       nfl=nf
 282       icg=mod(nf,2)+1
 283 cd      print *,'func',nf,nfl,icg
 284       call var_to_geom(n,x)
 285       call zerograd
 286       call chainbuild_extconf
 287 cd    write (iout,*) 'ETOTAL called from FUNC'
 288       call etotal(energia(0))
 289       call sum_gradient
 290       f=energia(0)
 291 c     if (jjj.gt.0) then
 292 c      write (iout,*) "upon exit from func"
 293 c      write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
 294 c      write (iout,*) 'f=',f
 295 c       jjj=0
 296 c     endif
 297       return
 298       end
 299 ************************************************************************
 300       subroutine func_restr(n,x,nf,f,uiparm,urparm,ufparm)
 301       implicit real*8 (a-h,o-z)
 302       include 'DIMENSIONS'
 303       include 'COMMON.DERIV'
 304       include 'COMMON.IOUNITS'
 305       include 'COMMON.GEO'
 306       common /chuju/ jjj
 307       double precision energia(0:n_ene)
 308       integer jjj
 309       double precision ufparm
 310       external ufparm
 311       integer uiparm(1)
 312       real*8 urparm(1)
 313       dimension x(maxvar)
 314 c     if (jjj.gt.0) then
 315 c       write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
 316 c     endif
 317       nfl=nf
 318       icg=mod(nf,2)+1
 319       call var_to_geom_restr(n,x)
 320       call zerograd
 321       call chainbuild_extconf
 322 cd    write (iout,*) 'ETOTAL called from FUNC'
 323       call etotal(energia(0))
 324       call sum_gradient
 325       f=energia(0)
 326 c     if (jjj.gt.0) then
 327 c       write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
 328 c       write (iout,*) 'f=',etot
 329 c       jjj=0
 330 c     endif
 331       return
 332       end
 333 c-------------------------------------------------------
 334       subroutine x2xx(x,xx,n)
 335       implicit real*8 (a-h,o-z)
 336       include 'DIMENSIONS'
 337       include 'COMMON.VAR'
 338       include 'COMMON.CHAIN'
 339       include 'COMMON.INTERACT'
 340       double precision xx(maxvar),x(maxvar)
 341
 342       do i=1,nvar
 343         varall(i)=x(i)
 344       enddo
 345
 346       ig=0
 347       igall=0
 348       do i=4,nres
 349         igall=igall+1
 350         if (mask_phi(i).eq.1) then
 351           ig=ig+1
 352           xx(ig)=x(igall)
 353         endif
 354       enddo
 355
 356       do i=3,nres
 357         igall=igall+1
 358         if (mask_theta(i).eq.1) then
 359           ig=ig+1
 360           xx(ig)=x(igall)
 361         endif
 362       enddo
 363
 364       do ij=1,2
 365       do i=2,nres-1
 366         if (itype(i).ne.10) then
 367           igall=igall+1
 368           if (mask_side(i).eq.1) then
 369             ig=ig+1
 370             xx(ig)=x(igall)
 371           endif
 372         endif
 373       enddo
 374       enddo
 375
 376       n=ig
 377
 378       return
 379       end
 380
 381       subroutine xx2x(x,xx)
 382       implicit real*8 (a-h,o-z)
 383       include 'DIMENSIONS'
 384       include 'COMMON.VAR'
 385       include 'COMMON.CHAIN'
 386       include 'COMMON.INTERACT'
 387       double precision xx(maxvar),x(maxvar)
 388
 389       do i=1,nvar
 390         x(i)=varall(i)
 391       enddo
 392
 393       ig=0
 394       igall=0
 395       do i=4,nres
 396         igall=igall+1
 397         if (mask_phi(i).eq.1) then
 398           ig=ig+1
 399           x(igall)=xx(ig)
 400         endif
 401       enddo
 402
 403       do i=3,nres
 404         igall=igall+1
 405         if (mask_theta(i).eq.1) then
 406           ig=ig+1
 407           x(igall)=xx(ig)
 408         endif
 409       enddo
 410
 411       do ij=1,2
 412       do i=2,nres-1
 413         if (itype(i).ne.10) then
 414           igall=igall+1
 415           if (mask_side(i).eq.1) then
 416             ig=ig+1
 417             x(igall)=xx(ig)
 418           endif
 419         endif
 420       enddo
 421       enddo
 422
 423       return
 424       end
 425
 426 c----------------------------------------------------------
 427       subroutine minim_dc(etot,iretcode,nfun)
 428 #ifdef LBFGS
 429       use minima
 430       use inform
 431       use output
 432       use iounit
 433       use scales
 434 #endif
 435       implicit real*8 (a-h,o-z)
 436       include 'DIMENSIONS'
 437 #ifndef LBFGS
 438       parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
 439 #endif
 440 #ifdef MPI
 441       include 'mpif.h'
 442 #endif
 443       include 'COMMON.CONTROL'
 444       include 'COMMON.SETUP'
 445       include 'COMMON.IOUNITS'
 446       include 'COMMON.VAR'
 447       include 'COMMON.GEO'
 448       include 'COMMON.MINIM'
 449       include 'COMMON.CHAIN'
 450       double precision minval,x(maxvar),d(maxvar),xx(maxvar)
 451 #ifdef LBFGS
 452       double precision grdmin
 453       double precision funcgrad_dc
 454       external funcgrad_dc
 455 #else
 456       dimension iv(liv)
 457       double precision v(1:lv)
 458       common /przechowalnia/ v
 459       external func_dc,grad_dc,fdum
 460 #endif
 461       double precision g(maxvar),f1
 462       integer nvarx
 463 #ifdef LBFGS
 464       maxiter=maxmin
 465       coordtype='CARTESIAN'
 466       grdmin=tolf
 467       jout=iout
 468       jprint=print_min_stat
 469       iwrite=0
 470 #else
 471       call deflt(2,iv,liv,lv,v)
 472 * 12 means fresh start, dont call deflt
 473       iv(1)=12
 474 * max num of fun calls
 475       if (maxfun.eq.0) maxfun=500
 476       iv(17)=maxfun
 477 * max num of iterations
 478       if (maxmin.eq.0) maxmin=1000
 479       iv(18)=maxmin
 480 * controls output
 481       iv(19)=2
 482 * selects output unit
 483       iv(21)=0
 484       if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout
 485 * 1 means to print out result
 486       iv(22)=print_min_res
 487 * 1 means to print out summary stats
 488       iv(23)=print_min_stat
 489 * 1 means to print initial x and d
 490       iv(24)=print_min_ini
 491 * min val for v(radfac) default is 0.1
 492       v(24)=0.1D0
 493 * max val for v(radfac) default is 4.0
 494       v(25)=2.0D0
 495 c     v(25)=4.0D0
 496 * check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
 497 * the sumsl default is 0.1
 498       v(26)=0.1D0
 499 * false conv if (act fnctn decrease) .lt. v(34)
 500 * the sumsl default is 100*machep
 501       v(34)=v(34)/100.0D0
 502 * absolute convergence
 503       if (tolf.eq.0.0D0) tolf=1.0D-4
 504       v(31)=tolf
 505 * relative convergence
 506       if (rtolf.eq.0.0D0) rtolf=1.0D-4
 507       v(32)=rtolf
 508 * controls initial step size
 509        v(35)=1.0D-1
 510 * large vals of d correspond to small components of step
 511       do i=1,6*nres
 512         d(i)=1.0D-1
 513       enddo
 514 #endif
 515       k=0
 516       do i=1,nres-1
 517         do j=1,3
 518           k=k+1
 519           x(k)=dc(j,i)
 520         enddo
 521       enddo
 522       do i=2,nres-1
 523         if (ialph(i,1).gt.0) then
 524         do j=1,3
 525           k=k+1
 526           x(k)=dc(j,i+nres)
 527         enddo
 528         endif
 529       enddo
 530       nvarx=k
 531       write (iout,*) "Variables set up nvarx",nvarx
 532 #ifdef LBFGS
 533 c
 534 c     From tinker
 535 c
 536 c     perform dynamic allocation of some global arrays
 537 c
 538       if (.not. allocated(scale))  allocate (scale(nvarx))
 539 c
 540 c     set scaling parameter for function and derivative values;
 541 c     use square root of median eigenvalue of typical Hessian
 542 c
 543       set_scale = .true.
 544 c      nvar = 0
 545       do i = 1, nvarx
 546 c         if (use(i)) then
 547 c            do j = 1, 3
 548 c               nvar = nvar + 1
 549                scale(i) = 12.0d0
 550 c            end do
 551 c         end if
 552       end do
 553       write (iout,*) "Calling lbfgs"
 554       call lbfgs (nvarx,x,etot,grdmin,funcgrad_dc,optsave)
 555       write (iout,*) "After lbfgs"
 556 #else
 557 c-----
 558 c      write (iout,*) "checkgrad before SUMSL"
 559 c      icheckgrad=1
 560 c      aincr=1.0d-7
 561 c      call exec_checkgrad
 562 c-----
 563
 564       call sumsl(k,d,x,func_dc,grad_dc,iv,liv,lv,v,idum,rdum,fdum)
 565 c-----
 566 c      write (iout,*) "checkgrad after SUMSL"
 567 c      call exec_checkgrad
 568 c-----
 569 #endif
 570       k=0
 571       do i=1,nres-1
 572         do j=1,3
 573           k=k+1
 574           dc(j,i)=x(k)
 575         enddo
 576       enddo
 577       do i=2,nres-1
 578         if (ialph(i,1).gt.0) then
 579         do j=1,3
 580           k=k+1
 581           dc(j,i+nres)=x(k)
 582         enddo
 583         endif
 584       enddo
 585       call chainbuild_cart
 586
 587 cd      call zerograd
 588 cd      nf=0
 589 cd      call func_dc(k,x,nf,f,idum,rdum,fdum)
 590 cd      call grad_dc(k,x,nf,g,idum,rdum,fdum)
 591 cd
 592 cd      do i=1,k
 593 cd       x(i)=x(i)+1.0D-5
 594 cd       call func_dc(k,x,nf,f1,idum,rdum,fdum)
 595 cd       x(i)=x(i)-1.0D-5
 596 cd       print '(i5,2f15.5)',i,g(i),(f1-f)/1.0D-5
 597 cd      enddo
 598 #ifndef LBFGS
 599       etot=v(10)
 600       iretcode=iv(1)
 601       nfun=iv(6)
 602 #endif
 603       return
 604       end
 605 #ifdef LBFGS
 606       double precision function funcgrad_dc(x,g)
 607       implicit real*8 (a-h,o-z)
 608       include 'DIMENSIONS'
 609 #ifdef MPI
 610       include 'mpif.h'
 611 #endif
 612       include 'COMMON.SETUP'
 613       include 'COMMON.CHAIN'
 614       include 'COMMON.DERIV'
 615       include 'COMMON.VAR'
 616       include 'COMMON.INTERACT'
 617       include 'COMMON.FFIELD'
 618       include 'COMMON.MD'
 619       include 'COMMON.IOUNITS'
 620       integer k
 621       dimension x(maxvar),g(maxvar)
 622       double precision energia(0:n_ene)
 623 c
 624       k=0
 625       do i=1,nres-1
 626         do j=1,3
 627           k=k+1
 628           dc(j,i)=x(k)
 629         enddo
 630       enddo
 631       do i=2,nres-1
 632         if (ialph(i,1).gt.0) then
 633         do j=1,3
 634           k=k+1
 635           dc(j,i+nres)=x(k)
 636         enddo
 637         endif
 638       enddo
 639       call chainbuild_cart
 640       call zerograd
 641       call etotal(energia(0))
 642 c      write (iout,*) "energia",energia(0)
 643       funcgrad_dc=energia(0)
 644 C
 645 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
 646 C
 647       call cartgrad
 648 cd      print *,40
 649       k=0
 650       do i=1,nres-1
 651         do j=1,3
 652           k=k+1
 653           g(k)=gcart(j,i)
 654         enddo
 655       enddo
 656       do i=2,nres-1
 657         if (ialph(i,1).gt.0) then
 658         do j=1,3
 659           k=k+1
 660           g(k)=gxcart(j,i)
 661         enddo
 662         endif
 663       enddo
 664
 665       return
 666       end
 667 #else
 668       subroutine func_dc(n,x,nf,f,uiparm,urparm,ufparm)
 669       implicit real*8 (a-h,o-z)
 670       include 'DIMENSIONS'
 671 #ifdef MPI
 672       include 'mpif.h'
 673 #endif
 674       include 'COMMON.SETUP'
 675       include 'COMMON.DERIV'
 676       include 'COMMON.IOUNITS'
 677       include 'COMMON.GEO'
 678       include 'COMMON.CHAIN'
 679       include 'COMMON.VAR'
 680       double precision energia(0:n_ene)
 681       double precision ufparm
 682       external ufparm
 683       integer uiparm(1)
 684       real*8 urparm(1)
 685       dimension x(maxvar)
 686       nfl=nf
 687 cbad      icg=mod(nf,2)+1
 688       icg=1
 689
 690       k=0
 691       do i=1,nres-1
 692         do j=1,3
 693           k=k+1
 694           dc(j,i)=x(k)
 695         enddo
 696       enddo
 697       do i=2,nres-1
 698         if (ialph(i,1).gt.0) then
 699         do j=1,3
 700           k=k+1
 701           dc(j,i+nres)=x(k)
 702         enddo
 703         endif
 704       enddo
 705       call chainbuild_cart
 706
 707       call zerograd
 708       call etotal(energia(0))
 709       f=energia(0)
 710
 711 cd      print *,'func_dc ',nf,nfl,f
 712
 713       return
 714       end
 715
 716       subroutine grad_dc(n,x,nf,g,uiparm,urparm,ufparm)
 717       implicit real*8 (a-h,o-z)
 718       include 'DIMENSIONS'
 719 #ifdef MPI
 720       include 'mpif.h'
 721 #endif
 722       include 'COMMON.SETUP'
 723       include 'COMMON.CHAIN'
 724       include 'COMMON.DERIV'
 725       include 'COMMON.VAR'
 726       include 'COMMON.INTERACT'
 727       include 'COMMON.FFIELD'
 728       include 'COMMON.MD'
 729       include 'COMMON.IOUNITS'
 730       external ufparm
 731       integer uiparm(1),k
 732       double precision urparm(1)
 733       dimension x(maxvar),g(maxvar)
 734 c
 735 c
 736 c
 737 cbad      icg=mod(nf,2)+1
 738       icg=1
 739 cd      print *,'grad_dc ',nf,nfl,nf-nfl+1,icg
 740       if (nf-nfl+1) 20,30,40
 741    20 call func_dc(n,x,nf,f,uiparm,urparm,ufparm)
 742 cd      print *,20
 743       if (nf.eq.0) return
 744       goto 40
 745    30 continue
 746 cd      print *,30
 747       k=0
 748       do i=1,nres-1
 749         do j=1,3
 750           k=k+1
 751           dc(j,i)=x(k)
 752         enddo
 753       enddo
 754       do i=2,nres-1
 755         if (ialph(i,1).gt.0) then
 756         do j=1,3
 757           k=k+1
 758           dc(j,i+nres)=x(k)
 759         enddo
 760         endif
 761       enddo
 762       call chainbuild_cart
 763
 764 C
 765 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
 766 C
 767    40 call cartgrad
 768 cd      print *,40
 769       k=0
 770       do i=1,nres-1
 771         do j=1,3
 772           k=k+1
 773           g(k)=gcart(j,i)
 774         enddo
 775       enddo
 776       do i=2,nres-1
 777         if (ialph(i,1).gt.0) then
 778         do j=1,3
 779           k=k+1
 780           g(k)=gxcart(j,i)
 781         enddo
 782         endif
 783       enddo
 784
 785       return
 786       end
 787 #endif