source/maxlik/src_MD_T_maxlik-NEWCORR-PMF-5/maxlikopt.F

   1       subroutine maxlikopt(nvarr,x)
   2       implicit none
   3       include "DIMENSIONS"
   4       include "DIMENSIONS.ZSCOPT"
   5 #ifdef MPI
   6       include "mpif.h"
   7       include "COMMON.MPI"
   8 #endif
   9       include "COMMON.NAMES"
  10       include "COMMON.WEIGHTS"
  11       include "COMMON.VMCPAR"
  12       include "COMMON.OPTIM"
  13       include "COMMON.CLASSES"
  14       include "COMMON.ENERGIES"
  15       include "COMMON.IOUNITS"
  16       integer i,j,k,iprot,nf,ii,inn,n,nn,indn(max_ene)
  17       logical solved,all_satisfied,not_done
  18       double precision ran_number,f,f0,x(max_paropt),x0(max_paropt),
  19      &  viol
  20       integer iran_num,nvarr
  21       double precision tcpu,t0,t1,t0w,t1w
  22       character*4 liczba
  23       integer ilen,lenpre
  24       external ilen
  25
  26       lenpre=ilen(prefix)
  27 #ifdef MPI
  28       if (me.eq.master) then
  29         if (nparmset.gt.0) then
  30           write (liczba,'(bz,i4.4)') myparm
  31           restartname=prefix(:lenpre)//"_par"//liczba//'.restart'
  32         else
  33           restartname=prefix(:lenpre)//'.restart'
  34         endif
  35       endif
  36 #else
  37       restartname=prefix(:lenpre)//'.restart'
  38 #endif
  39       print *,"Start solving inequalities"
  40       print *,"Mode",mode
  41       if (mode.eq.1) then
  42 c Only evaluate initial energies
  43         t0=tcpu()
  44 #ifdef MPI
  45         t0w=MPI_Wtime()
  46 #endif
  47         call minimize_setup(nvarr,x)
  48         t1=tcpu()
  49 #ifdef MPI
  50         write (iout,*) "CPU time for setup:",t1-t0," sec."
  51         t1w=MPI_Wtime()
  52         write (iout,*) "Wall clock time for setup:",t1w-t0w," sec."
  53         t0w=t1w
  54 #endif
  55         write (iout,*)
  56         do i=1,nvarr
  57           x0(i)=x(i)
  58         enddo
  59         call print_minimize_results(nvarr)
  60         return
  61       endif
  62       if (mode.ge.3) then
  63         t0=tcpu()
  64 #ifdef MPI
  65         t0w=MPI_Wtime()
  66 #endif
  67         call minimize_setup(nvarr,x)
  68         t1=tcpu()
  69 #ifdef MPI
  70         write (iout,*) "CPU time for setup:",t1-t0," sec."
  71         t1w=MPI_Wtime()
  72         write (iout,*) "Wall clock time for setup:",t1w-t0w," sec."
  73         t0w=t1w
  74 #endif
  75         t0=t1
  76         print *,"Calling minimize"
  77         call minimize(nvarr,x,f0)
  78         t1=tcpu()
  79         write (iout,*) "CPU time for minimization:",t1-t0," sec."
  80 #ifdef MPI
  81         t1w=MPI_Wtime()
  82         write (iout,*) "Wall clock time for minimization:",t1w-t0w
  83 #endif
  84         write (iout,*)
  85         do i=1,nvarr
  86           x0(i)=x(i)
  87         enddo
  88 c        if (nmcm.gt.0) then
  89 c          write (iout,*) "Start MCM procedure"
  90 c          call mcm
  91 c        endif
  92         call print_minimize_results(nvarr)
  93         return
  94       endif
  95 c
  96 c Scan weight space
  97 c
  98       if (mode.eq.2) then
  99         print *,"Calling scan"
 100         call minimize_setup(nvarr,x)
 101         call scan(nvarr,x,f0,.true.)
 102         do i=1,nvarr
 103           xm(i)=x(i)
 104         enddo
 105         call print_minimize_results(nvarr)
 106       endif
 107 50    format(bz,15f8.5)
 108 60    format(bz,100f10.5)
 109       return
 110       end
 111 c-----------------------------------------------------------------------------
 112       subroutine minimize_setup(nvar,x)
 113       implicit none
 114       include "DIMENSIONS"
 115       include "DIMENSIONS.ZSCOPT"
 116 #ifdef MPI
 117       include "mpif.h"
 118       integer IERROR,TAG,STATUS(MPI_STATUS_SIZE)
 119       include "COMMON.MPI"
 120 #endif
 121       include "COMMON.IOUNITS"
 122       include "COMMON.CLASSES"
 123       include "COMMON.NAMES"
 124       include "COMMON.OPTIM"
 125       include "COMMON.WEIGHTS"
 126       include "COMMON.WEIGHTDER"
 127       include "COMMON.VMCPAR"
 128       include "COMMON.ENERGIES"
 129       include "COMMON.TIME1"
 130       double precision viol
 131       integer i,j,k,ii,ib,ibatch,iprot,nvar,nf
 132       double precision x(max_paropt)
 133       logical lprn
 134       lprn=.true.
 135 #ifdef DEBUG
 136       write (iout,*) "Entered MINIMIZE_SETUP"
 137       call flush(iout)
 138 #endif
 139       call func1(nvar,3,x_orig)
 140       do i=1,n_ene
 141         ww_orig(i)=ww(i)
 142       enddo
 143       call x2w(nvar,x)
 144 c      write (iout,*) "Initial X",(x(i),i=1,nvar)
 145 c      write (iout,*) "X_orig",(x_orig(i),i=1,nvar)
 146 c      write (iout,*) "ww_oorig",(ww_oorig(i),i=1,n_ene)
 147 c      call flush(iout)
 148       call func1(nvar,3,x)
 149 #ifdef DEBUG
 150       write (iout,*) "x",(x(i),i=1,nvar)
 151 #endif
 152       write (iout,*)
 153       write(iout,'(10(1x,a6,1x))')(wname(print_order(i))(:6),i=1,n_ene)
 154       write(iout,40)(ww(print_order(i)),i=1,n_ene)
 155       write(iout,*)'-----------------------------------'
 156       call write_zscore
 157       call flush(iout)
 158       return
 159 40    format(10(f7.4,1x))
 160       end
 161 c------------------------------------------------------------------------------
 162       subroutine print_minimize_results(nvarr)
 163       implicit none
 164       include "DIMENSIONS"
 165       include "DIMENSIONS.ZSCOPT"
 166 #ifdef MPI
 167       include "mpif.h"
 168       integer ierror,status(MPI_STATUS_SIZE),tag
 169       include "COMMON.MPI"
 170 #endif
 171       include "COMMON.CONTROL"
 172       include "COMMON.WEIGHTS"
 173       include "COMMON.WEIGHTDER"
 174       include "COMMON.PROTNAME"
 175       include "COMMON.ENERGIES"
 176       include "COMMON.CLASSES"
 177       include "COMMON.IOUNITS"
 178       include "COMMON.VMCPAR"
 179       include "COMMON.OPTIM"
 180       include "COMMON.TIME1"
 181       integer nf,nvarr,i,ib,ic,ii,num_nat_tot,iprot,igather
 182       double precision x(max_ene)
 183       double precision viol
 184       character*4 liczba
 185       logical lprn
 186       lprn=.true.
 187 c      write (iout,*) "print_minimize_results"
 188       call flush(iout)
 189       write(iout,*) '----------------------'
 190       cutoffeval=.false.
 191       call x2w(nvarr,xm(1))
 192 c      write (iout,*) "Calling func1"
 193       call flush(iout)
 194       call func1(nvarr,2,xm(1))
 195 c      write (iout,*) "Calling write_params"
 196       call flush(iout)
 197       call write_params(nvarr,nf,xm(1))
 198 c      write(iout,*) "After write_params"
 199       call flush(iout)
 200       if (out_newe) then
 201       do ii=1,nprot
 202         if (nparmset.eq.1) then
 203         call write_prot(ii,'NewE_all')
 204         else
 205         write (liczba,'(bz,i4.4)') myparm
 206         call write_prot(ii,'NewE_all'//liczba)
 207         endif
 208 c        do ic=1,nclass(ii)
 209 c          call make_distrib(iclass(ic,ii),ii)
 210 c        enddo
 211       enddo
 212       endif
 213
 214       return
 215
 216       end
 217 c------------------------------------------------------------------------------
 218       subroutine write_params(nvar,nf,x)
 219       implicit none
 220       include "DIMENSIONS"
 221       include "DIMENSIONS.ZSCOPT"
 222       include "COMMON.CONTROL"
 223       include "COMMON.WEIGHTS"
 224       include "COMMON.NAMES"
 225       include "COMMON.INTERACT"
 226       include "COMMON.TORSION"
 227       include "COMMON.LOCAL"
 228       include "COMMON.SCCOR"
 229       include "COMMON.CLASSES"
 230       include "COMMON.ENERGIES"
 231       include "COMMON.IOUNITS"
 232       include "COMMON.FFIELD"
 233       include "COMMON.OPTIM"
 234       integer nf,nvar
 235       double precision x(max_paropt)
 236       integer ilen
 237       external ilen
 238       integer i,ii,j,l,ll,jj,pj,jstart,jend,it1,it2,k,iblock
 239       integer itor2typ(-2:2) /-20,9,10,9,20/
 240       character*1 aster(0:1) /" ","*"/,ast11,ast12,ast21,ast22
 241       character*4 liczba
 242       character*1 toronelet(-2:2) /"p","a","G","A","P"/
 243       logical lprint
 244       call x2w(nvar,x)
 245       call write_zscore
 246       write (iout,*) "Weights (asterisk: optimized)"
 247       write(iout,'(1x,15(a8,2x))') (wname(print_order(i))(:8),i=1,n_ene)
 248       write(iout,50)(ww(print_order(i)),
 249      &  aster(imask(print_order(i))),i=1,n_ene)
 250       write(iout,*)
 251 C Write weights in UNRES format for the input file
 252       if (nparmset.eq.1) then
 253         open(88,file="weights_opt."//prefix(:ilen(prefix)),
 254      &   status="unknown")
 255       else
 256         write (liczba,'(bz,i4.4)') myparm
 257         open(88,file="weights_opt."//prefix(:ilen(prefix))//"par_"
 258      &   //liczba,
 259      &   status="unknown")
 260       endif
 261 c      do i=1,n_ene,5
 262       do i=1,n_ene,4
 263         jstart=i
 264 c        jend=min0(i+4,n_ene)
 265         jend=min0(i+3,n_ene)
 266         jj=0
 267         do j=jstart,jend
 268           pj=print_order(j)
 269 c          write(88,'(bz,2a,f7.5," ",$)')
 270 c     &     wname(pj)(:ilen(wname(pj))),'=',ww(pj)
 271 c          jj=jj+ilen(wname(pj))+9
 272           write(88,'(bz,2a,e11.5," ",$)')
 273      &     wname(pj)(:ilen(wname(pj))),'=',ww(pj)
 274           jj=jj+ilen(wname(pj))+13
 275         enddo
 276         write (88,'(a,$)') (" ",j=1,79-jj)
 277         write (88,'("&")')
 278       enddo
 279       write (88,'(bz,2(2a,f7.5," "))')
 280      &    "CUTOFF","=",7.0d0,"WCORR4","=",0.0d0
 281       close(88)
 282       if (mod_side .or. mod_side_other) then
 283       if (nparmset.eq.1) then
 284       open(88,
 285      &     file="sc_"//POT(:ilen(POT))//"_opt."//prefix(:ilen(prefix)),
 286      &     status="unknown")
 287       else
 288       open(88,
 289      &     file="sc_"//POT(:ilen(POT))//"_opt."//prefix(:ilen(prefix))//
 290      &     "_par"//liczba,
 291      &     status="unknown")
 292       endif
 293       lprint=.true.
 294       write(88,'(2i5)') ipot,expon
 295       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
 296      & ', exponents are ',expon,2*expon
 297       goto (10,20,30,30,40) ipot
 298 C----------------------- LJ potential ---------------------------------
 299    10 write (88,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
 300       if (lprint) then
 301         write (iout,'(/a/)') 'Parameters of the LJ potential:'
 302         write (iout,'(a/)') 'The epsilon array:'
 303         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 304         write (iout,'(/a)') 'One-body parameters:'
 305         write (iout,'(a,4x,a)') 'residue','sigma'
 306         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
 307       endif
 308       goto 55
 309 C----------------------- LJK potential --------------------------------
 310    20 do i=1,ntyp
 311         write (88,'(4f20.10)')(eps(i,j),j=i,ntyp)
 312         write (iout,*)
 313       enddo
 314       write (88,'(4f20.10)') (sigma0(i),i=1,ntyp)
 315       write (88,*)
 316       write (88,'(4f20.10)') (rr0(i),i=1,ntyp)
 317       write (88,*)
 318       if (lprint) then
 319         write (iout,'(/a/)') 'Parameters of the LJK potential:'
 320         write (iout,'(a/)') 'The epsilon array:'
 321         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 322         write (iout,'(/a)') 'One-body parameters:'
 323         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
 324         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
 325      &        i=1,ntyp)
 326       endif
 327       goto 55
 328 C---------------------- GB or BP potential -----------------------------
 329    30 do i=1,ntyp
 330         write (88,'(4f20.15)') (eps(i,j),j=i,ntyp)
 331         write (88,*)
 332       enddo
 333       write (88,'(4f20.15)')(sigma0(i),i=1,ntyp)
 334       write (88,*)
 335       write (88,'(4f20.15)')(sigii(i),i=1,ntyp)
 336       write (88,*)
 337       write (88,'(4f20.15)')(chip0(i),i=1,ntyp)
 338       write (88,*)
 339       write (88,'(4f20.15)')(alp(i),i=1,ntyp)
 340       write (88,*)
 341 C For the GB potential convert sigma'**2 into chi'
 342       if (ipot.eq.4) then
 343         do i=1,ntyp
 344           chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
 345         enddo
 346       endif
 347       if (lprint) then
 348         if (ipot.eq.3) then
 349           write (iout,'(/a/)') 'Parameters of the BP potential:'
 350         else
 351           write (iout,'(/a/)') 'Parameters of the GB potential:'
 352         endif
 353         write (iout,'(a/)') 'The epsilon array:'
 354         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 355         write (iout,'(/a)') 'One-body parameters (asterisk: optimized):'
 356         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
 357      &       '    chip  ','    alph  '
 358         do i=1,ntyp
 359         write (iout,'(a3,6x,f10.5,1x,a1,3(f8.5,1x,a1))') restyp(i),
 360      & sigma0(i),aster(mask_sigma(1,i)),sigii(i),aster(mask_sigma(2,i)),
 361      &  chip(i),aster(mask_sigma(3,i)),alp(i),aster(mask_sigma(4,i))
 362         enddo
 363       endif
 364       goto 55
 365 C--------------------- GBV potential -----------------------------------
 366    40 write (88,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 367      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
 368      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
 369       if (lprint) then
 370         write (iout,'(/a/)') 'Parameters of the GBV potential:'
 371         write (iout,'(a/)') 'The epsilon array:'
 372         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 373         write (iout,'(/a)') 'One-body parameters:'
 374         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
 375      &      's||/s_|_^2','    chip  ','    alph  '
 376         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
 377      &           sigii(i),chip(i),alp(i),i=1,ntyp)
 378       endif
 379    55 continue
 380 C-----------------------------------------------------------------------
 381       close(88)
 382       endif
 383 #ifdef NEWCORR
 384       if (mod_tor .or. tor_mode.eq.2) then
 385 #else
 386       if (mod_tor) then
 387 #endif
 388         write (iout,'(a)') "Optimized weights of the torsional terms."
 389         do iblock=1,2
 390         do i=-ntortyp+1,ntortyp-1
 391           do j=-ntortyp+1,ntortyp-1
 392             if (mask_tor(0,j,i,iblock).gt.0) then
 393               write (iout,'(i4,2a4,f10.5)') 0,restyp(iloctyp(i)),
 394      &           restyp(iloctyp(j)),weitor(0,j,i,iblock)
 395             endif
 396           enddo
 397         enddo
 398         enddo
 399         if (nparmset.eq.1) then
 400
 401         open(88,file="tor_opt.parm."//prefix(:ilen(prefix)),
 402      &   status="unknown")
 403         write (88,'(i1,7h  ***  ,2a,7h  ***  )') ntortyp,
 404      &   "Optimized torsional parameters ",prefix(:ilen(prefix))
 405         write (88,'(40i2)') (itortyp(i),i=1,ntyp)
 406
 407         if (tor_mode.eq.0) then
 408
 409         write (iout,'(/a/)') 'Torsional constants:'
 410         do iblock=1,2
 411         do i=0,ntortyp-1
 412           do j=0,ntortyp-1
 413             write (iout,*) 'ityp',i,' jtyp',j
 414             write (iout,*) 'Fourier constants'
 415             do k=1,nterm(i,j,iblock)
 416               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
 417      &        v2(k,i,j,iblock)
 418             enddo
 419             write (iout,*) 'Lorenz constants'
 420             do k=1,nlor(i,j,iblock)
 421               write (iout,'(3(1pe15.5))')
 422      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
 423             enddo
 424           enddo
 425         enddo
 426         enddo
 427
 428         do iblock=1,2
 429         do i=0,ntortyp-1
 430           do j=-ntortyp+1,ntortyp-1
 431             write (88,'(2i2,16h  ************  ,a1,1h-,a1)')
 432      &       nterm(i,j,iblock),nlor(i,j,iblock),onelet(iloctyp(i)),
 433      &       onelet(iloctyp(j))
 434             do k=1,nterm(i,j,iblock)
 435               write (88,'(i6,13x,2(1pe18.5))') k,
 436      &         v1(k,i,j,iblock)*weitor(0,i,j,iblock),v2(k,i,j,
 437      &         iblock)*weitor(0,i,j,iblock)
 438             enddo
 439             do k=1,nlor(i,j,iblock)
 440               write (88,'(i6,13x,3(1pe18.5))')
 441      &         k,vlor1(k,i,j)*weitor(k,i,j,iblock),
 442      &         vlor2(k,i,j)*weitor(k,i,j,iblock),
 443      &         vlor3(k,i,j)*weitor(k,i,j,iblock)
 444             enddo
 445           enddo
 446         enddo
 447         enddo
 448
 449         else
 450 c Print valence-torsional parameters
 451         write (iout,'(a)')
 452      &    "Parameters of the valence-torsional potentials"
 453         do i=-ntortyp+1,ntortyp-1
 454         do j=-ntortyp+1,ntortyp-1
 455         write (iout,'(3a)')"Type ",onelet(iloctyp(i)),onelet(iloctyp(j))
 456         write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
 457         do k=1,nterm_kcc(j,i)
 458           do l=1,nterm_kcc_Tb(j,i)
 459             do ll=1,nterm_kcc_Tb(j,i)
 460                write (iout,'(3i5,2f15.4)')
 461      &            k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
 462             enddo
 463           enddo
 464         enddo
 465         enddo
 466         enddo
 467
 468         do i=-ntortyp+1,ntortyp-1
 469           do j=-ntortyp+1,ntortyp-1
 470             write (88,'(16h  ************  ,a1,1h-,a1)')
 471      &       onelet(iloctyp(i)),onelet(iloctyp(j))
 472             write (88,'(2i5)') nterm_kcc(j,i),nterm_kcc_Tb(j,i)
 473             do k=1,nterm_kcc(j,i)
 474               do l=1,nterm_kcc_Tb(j,i)
 475                 do ll=1,nterm_kcc_Tb(j,i)
 476                   write (88,'(3i5,2(1pe15.5))') k,l-1,ll-1,
 477      &            v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
 478                 enddo
 479               enddo
 480             enddo
 481           enddo
 482         enddo
 483
 484         close(88)
 485
 486         endif
 487
 488         endif
 489       endif
 490       if (mod_sccor) then
 491         write (iout,'(a)') "Optimized weights of the sccor terms."
 492         do i=1,nsccortyp
 493           do j=1,nsccortyp
 494             do l=1,3
 495             if (mask_tor(l,j,i,iblock).gt.0) then
 496               write (iout,'(i4,2a4,f10.5)') l,restyp(isccortyp(j)),
 497      &           restyp(isccortyp(i)),weitor(l,j,i,1)
 498             endif
 499             enddo
 500           enddo
 501         enddo
 502
 503         if (nparmset.eq.1) then
 504
 505         open(88,file="sccor_opt.parm."//prefix(:ilen(prefix)),
 506      &     status="unknown")
 507         write (88,'(i2,7h  ***  ,2a,7h  ***  )') nsccortyp,
 508      & "Optimized local correlation parameters ",prefix(:ilen(prefix))
 509         write (88,'(20i4)') (isccortyp(i),i=1,ntyp)
 510         do l=1,3
 511         do i=1,nsccortyp
 512           do j=1,nsccortyp
 513             write (88,'(2i4)') nterm_sccor(i,j),nlor_sccor(i,j)
 514             write (88,'(12(1h*),1x,a1,1h-,a1)') onelet(i),onelet(j)
 515             do k=1,nterm_sccor(i,j)
 516               write (88,'(i6,13x,2(1pe18.5))') k,
 517      &          v1sccor(k,l,i,j)*weitor(l,i,j,1),
 518      &          v2sccor(k,l,i,j)*weitor(l,i,j,1)
 519             enddo
 520           enddo
 521         enddo
 522         enddo
 523         close(88)
 524
 525         endif
 526
 527       endif
 528 C-----------------------------------------------------------------------
 529 c 7/8/17 AL: Optimization of the bending parameters
 530       if (mod_ang) then
 531         write(iout,*)
 532      &    "Optimized angle potential coefficients (same for D-types)"
 533         do i=0,nthetyp-1
 534           if (mask_ang(i).eq.0) cycle
 535           write (iout,*) "Type: ",onelet(iloctyp(i))
 536           do j=1,nbend_kcc_TB(i)
 537             write (iout,'(i5,f15.5)') j,v1bend_chyb(j,i)
 538           enddo
 539         enddo
 540         open(88,file="theta_opt.parm."//prefix(:ilen(prefix)),
 541      &   status="unknown")
 542         write (88,'(i2)') nthetyp
 543         do i=-nthetyp+1,nthetyp-1
 544           write (88,'(i2,1x,12(1h*),1x,a1)')
 545      &      nbend_kcc_TB(i),onelet(iloctyp(i))
 546           do j=0,nbend_kcc_TB(i)
 547             write (88,'(i5,f20.10)') j,v1bend_chyb(j,i)
 548           enddo
 549         enddo
 550         close(88)
 551       endif
 552 C-----------------------------------------------------------------------
 553       if (mod_fourier(nloctyp).gt.0) then
 554 #ifdef NEWCORR
 555         write (iout,'(2a)')
 556      &   "Fourier coefficients of new (angle-dependent) cumulants",
 557      &   " (asterisk: optimized)"
 558         do i=0,nloctyp-1
 559           write (iout,*) "Type: ",onelet(iloctyp(i))
 560           write (iout,*) "Coefficients of the expansion of B1"
 561           do j=1,2
 562             write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
 563           enddo
 564           write (iout,*) "Coefficients of the expansion of B2"
 565           do j=1,2
 566             write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
 567           enddo
 568           write (iout,*) "Coefficients of the expansion of C"
 569           write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
 570           write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
 571           write (iout,*) "Coefficients of the expansion of D"
 572           write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
 573           write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
 574           write (iout,*) "Coefficients of the expansion of E"
 575           write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
 576           do j=1,2
 577             do k=1,2
 578               write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
 579             enddo
 580           enddo
 581         enddo
 582
 583         IF (SPLIT_FOURIERTOR) THEN
 584
 585         write (iout,'(2a)')
 586      &   "Fourier coefficients of new (angle-dependent) cumulants",
 587      &   " for torsional potentials (asterisk: optimized)"
 588         do i=0,nloctyp-1
 589           write (iout,*) "Type: ",onelet(iloctyp(i))
 590           write (iout,*) "Coefficients of the expansion of B1"
 591           do j=1,2
 592             write (iout,'(3hB1(,i1,1h),3f10.5)')
 593      &         j,(bnew1tor(k,j,i),k=1,3)
 594           enddo
 595           write (iout,*) "Coefficients of the expansion of B2"
 596           do j=1,2
 597             write (iout,'(3hB2(,i1,1h),3f10.5)')
 598      &         j,(bnew2tor(k,j,i),k=1,3)
 599           enddo
 600           write (iout,*) "Coefficients of the expansion of C"
 601           write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
 602           write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
 603           write (iout,*) "Coefficients of the expansion of D"
 604           write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
 605           write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
 606           write (iout,*) "Coefficients of the expansion of E"
 607           write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
 608           do j=1,2
 609             do k=1,2
 610               write (iout,'(1hE,2i1,2f10.5)')
 611      &         j,k,(eenewtor(l,j,k,i),l=1,2)
 612             enddo
 613           enddo
 614         enddo
 615
 616         ENDIF
 617
 618 #else
 619         write (iout,'(2a)')
 620      &   "Fourier coefficients of old angle-independent cumulants",
 621      &   " (asterisk: optimized)"
 622         do i=0,nloctyp-1
 623           write (iout,*) 'Type ',restyp(iloctyp(i))
 624             write (iout,'(a,i2,a,f10.5,1x,a)')
 625      &    ('b(',j,')=',b(j,i),aster(mask_fourier(j,i)),j=1,13)
 626         enddo
 627 #endif
 628 c Write a file with parameters for UNRES
 629         if (nparmset.eq.1) then
 630         open(88,file="fourier_opt.parm."//prefix(:ilen(prefix)),
 631      &     status="unknown")
 632         else
 633         open(88,file="fourier_opt.parm."//prefix(:ilen(prefix))//
 634      &     "_par"//liczba,
 635      &     status="unknown")
 636         endif
 637 #ifdef NEWCORR
 638         IF (SPLIT_FOURIERTOR) THEN
 639           write (88,'(i5,5x,a)') -nloctyp,
 640      &                         "# Number of local interaction types"
 641         ELSE
 642           write (88,'(i5,5x,a)') nloctyp,
 643      &                         "# Number of local interaction types"
 644         ENDIF
 645         write (88,'(40i2)') (itype2loc(i),i=1,ntyp)
 646         write (88,'(20i4)') (iloctyp(i),i=0,nloctyp)
 647         do i=0,nloctyp-1
 648           write (88,'(a)') restyp(iloctyp(i))
 649           do ii=1,3
 650             do j=1,2
 651               write (88,'(f20.10,4h  b1,i1,1h(,i1,1h))')
 652      &         bnew1(ii,j,i),j,ii
 653             enddo
 654           enddo
 655           do ii=1,3
 656             do j=1,2
 657               write (88,'(f20.10,4h  b2,i1,1h(,i1,1h))')
 658      &         bnew2(ii,j,i),j,ii
 659             enddo
 660           enddo
 661           do ii=1,3
 662             do j=1,2
 663               write (88,'(f20.10,4h  c1,i1,1h(,i1,1h))')
 664      &         2*ccnew(ii,j,i),j,ii
 665             enddo
 666           enddo
 667           do ii=1,3
 668             do j=1,2
 669               write (88,'(f20.10,4h  d1,i1,1h(,i1,1h))')
 670      &         2*ddnew(ii,j,i),j,ii
 671             enddo
 672           enddo
 673           do ii=1,2
 674             do j=1,2
 675               do k=1,2
 676               write (88,'(f20.10,3h  e,2i1,1h(,i1,1h))')
 677      &         eenew(ii,j,k,i),j,k,ii
 678               enddo
 679             enddo
 680           enddo
 681           do ii=1,3
 682               write (88,'(f20.10,5h  e0(,i1,1h))')
 683      &         e0new(ii,i),ii
 684           enddo
 685         enddo
 686
 687         IF (SPLIT_FOURIERTOR) THEN
 688
 689         do i=0,nloctyp-1
 690           write (88,'(a)') restyp(iloctyp(i))
 691           do ii=1,3
 692             do j=1,2
 693               write (88,'(f20.10,4h  b1,i1,1h(,i1,1h))')
 694      &         bnew1tor(ii,j,i),j,ii
 695             enddo
 696           enddo
 697           do ii=1,3
 698             do j=1,2
 699               write (88,'(f20.10,4h  b2,i1,1h(,i1,1h))')
 700      &         bnew2tor(ii,j,i),j,ii
 701             enddo
 702           enddo
 703           do ii=1,3
 704             do j=1,2
 705               write (88,'(f20.10,4h  c1,i1,1h(,i1,1h))')
 706      &         2*ccnewtor(ii,j,i),j,ii
 707             enddo
 708           enddo
 709           do ii=1,3
 710             do j=1,2
 711               write (88,'(f20.10,4h  d1,i1,1h(,i1,1h))')
 712      &         2*ddnewtor(ii,j,i),j,ii
 713             enddo
 714           enddo
 715           do ii=1,2
 716             do j=1,2
 717               do k=1,2
 718               write (88,'(f20.10,3h  e,2i1,1h(,i1,1h))')
 719      &         eenewtor(ii,j,k,i),j,k,ii
 720               enddo
 721             enddo
 722           enddo
 723           do ii=1,3
 724               write (88,'(f20.10,5h  e0(,i1,1h))')
 725      &         e0newtor(ii,i),ii
 726           enddo
 727         enddo
 728
 729         ENDIF
 730
 731 #else
 732         write (88,'(i5,5x,a)') nloctyp,
 733      &                         "# Number of local interaction types"
 734         do i=1,nloctyp
 735           write (88,'(a)') restyp(iloctyp(i))
 736           do j=1,13
 737             write (88,'(f20.15)') b(j,i)
 738           enddo
 739         enddo
 740 #endif
 741         close(88)
 742       endif
 743       if (mod_elec) then
 744                 write (iout,'(/a)')
 745      &  'Electrostatic interaction constants (asterisk: optimized):'
 746         write (iout,'(1x,a,1x,a,7x,a,15x,a,15x,a,13x,a)')
 747      &            'IT','JT','EPP','RPP','ELPP6','ELPP3'
 748         do i=1,2
 749           do j=1,2
 750             write(iout,'(2i3,4(1pe15.4,1x,a1,1x))')i,j,
 751      &       epp(i,j),aster(mask_elec(i,j,1)),
 752      &       rpp(i,j),aster(mask_elec(i,j,2)),
 753      &       elpp6(i,j),aster(mask_elec(i,j,3)),
 754      &       elpp3(i,j),aster(mask_elec(i,j,4))
 755           enddo
 756         enddo
 757         write (iout,*)
 758         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
 759      &            'IT','JT','APP','BPP','AEL6','AEL3'
 760         do i=1,2
 761           do j=1,2
 762             write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
 763      &                    ael6(i,j),ael3(i,j)
 764           enddo
 765         enddo
 766 C Write electrostatic interaction parameters to a file
 767         if (nparmset.eq.1) then
 768         open(88,file="electr_opt.parm."//prefix(:ilen(prefix)),
 769      &   status="unknown")
 770         else
 771         open(88,file="electr_opt.parm."//prefix(:ilen(prefix))//
 772      &   "_par"//liczba,
 773      &   status="unknown")
 774         endif
 775         write(88,'(4f15.10,5x,a)')((epp(i,j),j=1,2),i=1,2),"! EPP"
 776         write(88,'(4f15.10,5x,a)')((rpp(i,j),j=1,2),i=1,2),"! RPP"
 777         write(88,'(4f15.10,5x,a)')((elpp6(i,j),j=1,2),i=1,2),"! ELPP6"
 778         write(88,'(4f15.10,5x,a)')((elpp3(i,j),j=1,2),i=1,2),"! ELPP3"
 779         close(88)
 780       endif
 781       if (mod_scp) then
 782         write (iout,*)
 783         write (iout,*)
 784      &  "Parameters of SC-p interactions (star: optimized):"
 785         write (iout,'(t14,a,t34,a)') "TYPE 1","TYPE 2"
 786         write (iout,'(8a)') "SC TYPE"," EPS_SCP  ","   R_SCP  ",
 787      &    " EPS_SCP  ","   R_SCP  "
 788         do i=1,20
 789           ast11 = aster(mask_scp(i,1,1))
 790           ast12 = aster(mask_scp(i,1,2))
 791           ast21 = aster(mask_scp(i,2,1))
 792           ast22 = aster(mask_scp(i,2,2))
 793           if (mask_scp(i,1,1).eq.0) ast11 = aster(mask_scp(0,1,1))
 794           if (mask_scp(i,1,2).eq.0) ast12 = aster(mask_scp(0,1,2))
 795           if (mask_scp(i,2,1).eq.0) ast21 = aster(mask_scp(0,2,1))
 796           if (mask_scp(i,2,2).eq.0) ast22 = aster(mask_scp(0,2,2))
 797           write (iout,'(2x,A3,2x,4(f8.3,1x,a1))') restyp(i),
 798      &     eps_scp(i,1),ast11,rscp(i,1),ast12,eps_scp(i,2),ast21,
 799      &     rscp(i,2),ast22
 800         enddo
 801 C Write SCp parameters to a file
 802         if (nparmset.eq.1) then
 803         open(88,file="scp_opt.parm."//prefix(:ilen(prefix)),
 804      &     status="unknown")
 805         else
 806         open(88,file="scp_opt.parm."//prefix(:ilen(prefix))//
 807      &     "_par"//liczba,
 808      &     status="unknown")
 809         endif
 810         do i=1,ntyp
 811           write(88,'(4f15.10,5x,a)') eps_scp(i,1),rscp(i,1),
 812      &      eps_scp(i,2),rscp(i,2),restyp(i)
 813         enddo
 814         close(88)
 815       endif
 816 50    format(bz,15(f8.5,1x,a1))
 817 60    format(bz,100f10.5)
 818       end
 819 c------------------------------------------------------------------------------
 820       subroutine write_zscore
 821       implicit none
 822       include "DIMENSIONS"
 823       include "DIMENSIONS.ZSCOPT"
 824       include "COMMON.ENERGIES"
 825       include "COMMON.CLASSES"
 826       include "COMMON.PROTNAME"
 827       include "COMMON.VMCPAR"
 828       include "COMMON.IOUNITS"
 829       include "COMMON.WEIGHTS"
 830       include "COMMON.WEIGHTDER"
 831       include "COMMON.COMPAR"
 832       include "COMMON.OPTIM"
 833       include "COMMON.THERMAL"
 834       integer i,j,k,iprot,ib
 835       character*6 namnat(5) /"angrms","Q","rmsd","rgy","sign"/
 836       integer ilen
 837       external ilen
 838       do iprot=1,nprot
 839
 840         write (iout,*)
 841         write (iout,'(a,2x,a)') "Protein",
 842      &    protname(iprot)(:ilen(protname(iprot)))
 843
 844         write (iout,*)
 845         write (iout,'(a)')"Maximum likelihood."
 846         write (iout,*)
 847         do ib=1,nbeta(1,iprot)
 848           write (iout,'(f8.1,f10.5)')
 849      &      1.0d0/(Rgas*betaT(ib,1,iprot)),sumlik(ib,iprot)
 850         enddo
 851         write (iout,*)
 852
 853         write (iout,*)
 854         write (iout,'(a)')"Specific heat."
 855         write (iout,*)
 856         write (iout,'(a,f10.5)') "Baseline",heatbase(iprot)
 857         do ib=1,nbeta(2,iprot)
 858           write (iout,'(f8.1,3f10.5)')
 859      &      1.0d0/(Rgas*betaT(ib,2,iprot)),
 860      &      zvtot(ib,iprot),target_cv(ib,iprot),weiCv(ib,iprot)
 861         enddo
 862
 863         if (natlike(iprot).gt.0) write (iout,'(a)')"Nativelikeness."
 864         do i=1,natlike(iprot)
 865           write (iout,*) "Property",i,numnat(i,iprot)
 866           if (natdim(i,iprot).eq.1)  then
 867             do ib=1,nbeta(i+2,iprot)
 868               write (iout,'(f7.2,2f10.5)')
 869      &         1.0d0/(Rgas*betaT(ib,i+2,iprot)),nuave(1,ib,i,iprot),
 870      &         nuexp(1,ib,i,iprot),weinat(1,ib,i,iprot)
 871             enddo
 872           else
 873             do ib=1,nbeta(i+2,iprot)
 874               write (iout,'("Temperature",f7.2," NATDIM",i5)')
 875      &        1.0d0/(Rgas*betaT(ib,i+2,iprot)),natdim(i,iprot)
 876               do k=1,natdim(i,iprot)
 877               write (iout,'(3f10.5)') nuave(k,ib,i,iprot),
 878      &         nuexp(k,ib,i,iprot),weinat(k,ib,i,iprot)
 879               enddo
 880             enddo
 881           endif
 882         enddo
 883
 884       enddo
 885       write (iout,*)
 886       return
 887       end
 888 c------------------------------------------------------------------------------
 889       subroutine write_prot(ii,przed)
 890       implicit none
 891       include "DIMENSIONS.ZSCOPT"
 892       include "COMMON.ENERGIES"
 893       include "COMMON.CLASSES"
 894       include "COMMON.PROTNAME"
 895       include "COMMON.VMCPAR"
 896       include "COMMON.IOUNITS"
 897       include "COMMON.COMPAR"
 898       integer i,ii,j
 899       integer ilen
 900       external ilen
 901       character*(*) przed
 902       character*64 nazwa
 903       nazwa=przed(:ilen(przed))//'.'//prefix(:ilen(prefix))
 904      &   //'.'//protname(ii)
 905       open(unit=istat,file=nazwa)
 906       do i=1,ntot_work(ii)
 907         write(istat,'(i10,3f15.3,20e15.5)')
 908      &    i,e_total(i,ii),eini(i,ii),entfac(i,ii),
 909      &    (Ptab(i,j,ii),j=1,nbeta(1,ii))
 910       enddo
 911       close(istat)
 912       return
 913       end