source/maxlik/src_MD_T_maxlik-NEWCORR.safe/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.SCCOR"
 228       include "COMMON.CLASSES"
 229       include "COMMON.ENERGIES"
 230       include "COMMON.IOUNITS"
 231       include "COMMON.FFIELD"
 232       include "COMMON.OPTIM"
 233       integer nf,nvar
 234       double precision x(max_paropt)
 235       integer ilen
 236       external ilen
 237       integer i,ii,j,l,ll,jj,pj,jstart,jend,it1,it2,k,iblock
 238       integer itor2typ(-2:2) /-20,9,10,9,20/
 239       character*1 aster(0:1) /" ","*"/,ast11,ast12,ast21,ast22
 240       character*4 liczba
 241       character*1 toronelet(-2:2) /"p","a","G","A","P"/
 242       logical lprint
 243       call x2w(nvar,x)
 244       call write_zscore
 245       write (iout,*) "Weights (asterisk: optimized)"
 246       write(iout,'(1x,15(a8,2x))') (wname(print_order(i))(:8),i=1,n_ene)
 247       write(iout,50)(ww(print_order(i)),
 248      &  aster(imask(print_order(i))),i=1,n_ene)
 249       write(iout,*)
 250 C Write weights in UNRES format for the input file
 251       if (nparmset.eq.1) then
 252         open(88,file="weights_opt."//prefix(:ilen(prefix)),
 253      &   status="unknown")
 254       else
 255         write (liczba,'(bz,i4.4)') myparm
 256         open(88,file="weights_opt."//prefix(:ilen(prefix))//"par_"
 257      &   //liczba,
 258      &   status="unknown")
 259       endif
 260 c      do i=1,n_ene,5
 261       do i=1,n_ene,4
 262         jstart=i
 263 c        jend=min0(i+4,n_ene)
 264         jend=min0(i+3,n_ene)
 265         jj=0
 266         do j=jstart,jend
 267           pj=print_order(j)
 268 c          write(88,'(bz,2a,f7.5," ",$)')
 269 c     &     wname(pj)(:ilen(wname(pj))),'=',ww(pj)
 270 c          jj=jj+ilen(wname(pj))+9
 271           write(88,'(bz,2a,e11.5," ",$)')
 272      &     wname(pj)(:ilen(wname(pj))),'=',ww(pj)
 273           jj=jj+ilen(wname(pj))+13
 274         enddo
 275         write (88,'(a,$)') (" ",j=1,79-jj)
 276         write (88,'("&")')
 277       enddo
 278       write (88,'(bz,2(2a,f7.5," "))')
 279      &    "CUTOFF","=",7.0d0,"WCORR4","=",0.0d0
 280       close(88)
 281       if (mod_side .or. mod_side_other) then
 282       if (nparmset.eq.1) then
 283       open(88,
 284      &     file="sc_"//POT(:ilen(POT))//"_opt."//prefix(:ilen(prefix)),
 285      &     status="unknown")
 286       else
 287       open(88,
 288      &     file="sc_"//POT(:ilen(POT))//"_opt."//prefix(:ilen(prefix))//
 289      &     "_par"//liczba,
 290      &     status="unknown")
 291       endif
 292       lprint=.true.
 293       write(88,'(2i5)') ipot,expon
 294       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
 295      & ', exponents are ',expon,2*expon
 296       goto (10,20,30,30,40) ipot
 297 C----------------------- LJ potential ---------------------------------
 298    10 write (88,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
 299       if (lprint) then
 300         write (iout,'(/a/)') 'Parameters of the LJ potential:'
 301         write (iout,'(a/)') 'The epsilon array:'
 302         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 303         write (iout,'(/a)') 'One-body parameters:'
 304         write (iout,'(a,4x,a)') 'residue','sigma'
 305         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
 306       endif
 307       goto 55
 308 C----------------------- LJK potential --------------------------------
 309    20 write (88,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 310      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
 311       if (lprint) then
 312         write (iout,'(/a/)') 'Parameters of the LJK potential:'
 313         write (iout,'(a/)') 'The epsilon array:'
 314         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 315         write (iout,'(/a)') 'One-body parameters:'
 316         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
 317         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
 318      &        i=1,ntyp)
 319       endif
 320       goto 55
 321 C---------------------- GB or BP potential -----------------------------
 322    30 do i=1,ntyp
 323         write (88,'(4f20.15)') (eps(i,j),j=i,ntyp)
 324       enddo
 325       write (88,'(4f20.15)')(sigma0(i),i=1,ntyp)
 326       write (88,'(4f20.15)')(sigii(i),i=1,ntyp)
 327       write (88,'(4f20.15)')(chip0(i),i=1,ntyp)
 328       write (88,'(4f20.15)')(alp(i),i=1,ntyp)
 329 C For the GB potential convert sigma'**2 into chi'
 330       if (ipot.eq.4) then
 331         do i=1,ntyp
 332           chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
 333         enddo
 334       endif
 335       if (lprint) then
 336         if (ipot.eq.3) then
 337           write (iout,'(/a/)') 'Parameters of the BP potential:'
 338         else
 339           write (iout,'(/a/)') 'Parameters of the GB potential:'
 340         endif
 341         write (iout,'(a/)') 'The epsilon array:'
 342         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 343         write (iout,'(/a)') 'One-body parameters (asterisk: optimized):'
 344         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
 345      &       '    chip  ','    alph  '
 346         do i=1,ntyp
 347         write (iout,'(a3,6x,f10.5,1x,a1,3(f8.5,1x,a1))') restyp(i),
 348      & sigma0(i),aster(mask_sigma(1,i)),sigii(i),aster(mask_sigma(2,i)),
 349      &  chip(i),aster(mask_sigma(3,i)),alp(i),aster(mask_sigma(4,i))
 350         enddo
 351       endif
 352       goto 55
 353 C--------------------- GBV potential -----------------------------------
 354    40 write (88,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 355      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
 356      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
 357       if (lprint) then
 358         write (iout,'(/a/)') 'Parameters of the GBV potential:'
 359         write (iout,'(a/)') 'The epsilon array:'
 360         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 361         write (iout,'(/a)') 'One-body parameters:'
 362         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
 363      &      's||/s_|_^2','    chip  ','    alph  '
 364         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
 365      &           sigii(i),chip(i),alp(i),i=1,ntyp)
 366       endif
 367    55 continue
 368 C-----------------------------------------------------------------------
 369       close(88)
 370       endif
 371 #ifdef NEWCORR
 372       if (mod_tor .or. tor_mode.eq.2) then
 373 #else
 374       if (mod_tor) then
 375 #endif
 376         write (iout,'(a)') "Optimized weights of the torsional terms."
 377         do iblock=1,2
 378         do i=-ntortyp+1,ntortyp-1
 379           do j=-ntortyp+1,ntortyp-1
 380             if (mask_tor(0,j,i,iblock).gt.0) then
 381               write (iout,'(i4,2a4,f10.5)') 0,restyp(iloctyp(i)),
 382      &           restyp(iloctyp(j)),weitor(0,j,i,iblock)
 383             endif
 384           enddo
 385         enddo
 386         enddo
 387         if (nparmset.eq.1) then
 388
 389         open(88,file="tor_opt.parm."//prefix(:ilen(prefix)),
 390      &   status="unknown")
 391         write (88,'(i1,7h  ***  ,2a,7h  ***  )') ntortyp,
 392      &   "Optimized torsional parameters ",prefix(:ilen(prefix))
 393         write (88,'(40i2)') (itortyp(i),i=1,ntyp)
 394
 395         if (tor_mode.eq.0) then
 396
 397         write (iout,'(/a/)') 'Torsional constants:'
 398         do iblock=1,2
 399         do i=0,ntortyp-1
 400           do j=0,ntortyp-1
 401             write (iout,*) 'ityp',i,' jtyp',j
 402             write (iout,*) 'Fourier constants'
 403             do k=1,nterm(i,j,iblock)
 404               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
 405      &        v2(k,i,j,iblock)
 406             enddo
 407             write (iout,*) 'Lorenz constants'
 408             do k=1,nlor(i,j,iblock)
 409               write (iout,'(3(1pe15.5))')
 410      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
 411             enddo
 412           enddo
 413         enddo
 414         enddo
 415
 416         do iblock=1,2
 417         do i=0,ntortyp-1
 418           do j=-ntortyp+1,ntortyp-1
 419             write (88,'(2i2,16h  ************  ,a1,1h-,a1)')
 420      &       nterm(i,j,iblock),nlor(i,j,iblock),onelet(iloctyp(i)),
 421      &       onelet(iloctyp(j))
 422             do k=1,nterm(i,j,iblock)
 423               write (88,'(i6,13x,2(1pe18.5))') k,
 424      &         v1(k,i,j,iblock)*weitor(0,i,j,iblock),v2(k,i,j,
 425      &         iblock)*weitor(0,i,j,iblock)
 426             enddo
 427             do k=1,nlor(i,j,iblock)
 428               write (88,'(i6,13x,3(1pe18.5))')
 429      &         k,vlor1(k,i,j)*weitor(k,i,j,iblock),
 430      &         vlor2(k,i,j)*weitor(k,i,j,iblock),
 431      &         vlor3(k,i,j)*weitor(k,i,j,iblock)
 432             enddo
 433           enddo
 434         enddo
 435         enddo
 436
 437         else
 438 c Print valence-torsional parameters
 439         write (iout,'(a)')
 440      &    "Parameters of the valence-torsional potentials"
 441         do i=-ntortyp+1,ntortyp-1
 442         do j=-ntortyp+1,ntortyp-1
 443         write (iout,'(3a)')"Type ",onelet(iloctyp(i)),onelet(iloctyp(j))
 444         write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
 445         do k=1,nterm_kcc(j,i)
 446           do l=1,nterm_kcc_Tb(j,i)
 447             do ll=1,nterm_kcc_Tb(j,i)
 448                write (iout,'(3i5,2f15.4)')
 449      &            k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
 450             enddo
 451           enddo
 452         enddo
 453         enddo
 454         enddo
 455
 456         do i=-ntortyp+1,ntortyp-1
 457           do j=-ntortyp+1,ntortyp-1
 458             write (88,'(16h  ************  ,a1,1h-,a1)')
 459      &       onelet(iloctyp(i)),onelet(iloctyp(j))
 460             write (88,'(2i5)') nterm_kcc(j,i),nterm_kcc_Tb(j,i)
 461             do k=1,nterm_kcc(j,i)
 462               do l=1,nterm_kcc_Tb(j,i)
 463                 do ll=1,nterm_kcc_Tb(j,i)
 464                   write (88,'(3i5,2(1pe15.5))') k,l-1,ll-1,
 465      &            v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
 466                 enddo
 467               enddo
 468             enddo
 469           enddo
 470         enddo
 471
 472         close(88)
 473
 474         endif
 475
 476         endif
 477       endif
 478       if (mod_sccor) then
 479         write (iout,'(a)') "Optimized weights of the sccor terms."
 480         do i=1,nsccortyp
 481           do j=1,nsccortyp
 482             do l=1,3
 483             if (mask_tor(l,j,i,iblock).gt.0) then
 484               write (iout,'(i4,2a4,f10.5)') l,restyp(isccortyp(j)),
 485      &           restyp(isccortyp(i)),weitor(l,j,i,1)
 486             endif
 487             enddo
 488           enddo
 489         enddo
 490
 491         if (nparmset.eq.1) then
 492
 493         open(88,file="sccor_opt.parm."//prefix(:ilen(prefix)),
 494      &     status="unknown")
 495         write (88,'(i2,7h  ***  ,2a,7h  ***  )') nsccortyp,
 496      & "Optimized local correlation parameters ",prefix(:ilen(prefix))
 497         write (88,'(20i4)') (isccortyp(i),i=1,ntyp)
 498         do l=1,3
 499         do i=1,nsccortyp
 500           do j=1,nsccortyp
 501             write (88,'(2i4)') nterm_sccor(i,j),nlor_sccor(i,j)
 502             write (88,'(12(1h*),1x,a1,1h-,a1)') onelet(i),onelet(j)
 503             do k=1,nterm_sccor(i,j)
 504               write (88,'(i6,13x,2(1pe18.5))') k,
 505      &          v1sccor(k,l,i,j)*weitor(l,i,j,1),
 506      &          v2sccor(k,l,i,j)*weitor(l,i,j,1)
 507             enddo
 508           enddo
 509         enddo
 510         enddo
 511         close(88)
 512
 513         endif
 514
 515       endif
 516 C-----------------------------------------------------------------------
 517       if (mod_fourier(nloctyp).gt.0) then
 518         write (iout,*)
 519      &   "Fourier coefficients of cumulants (asterisk: optimized)"
 520         do i=0,nloctyp-1
 521 #ifdef NEWCORR
 522           write (iout,*) "Type: ",onelet(iloctyp(i))
 523           write (iout,*) "Coefficients of the expansion of B1"
 524           do j=1,2
 525             write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
 526           enddo
 527           write (iout,*) "Coefficients of the expansion of B2"
 528           do j=1,2
 529             write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
 530           enddo
 531           write (iout,*) "Coefficients of the expansion of C"
 532           write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
 533           write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
 534           write (iout,*) "Coefficients of the expansion of D"
 535           write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
 536           write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
 537           write (iout,*) "Coefficients of the expansion of E"
 538           write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
 539           do j=1,2
 540             do k=1,2
 541               write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
 542             enddo
 543           enddo
 544 #else
 545           write (iout,*) 'Type ',restyp(iloctyp(i))
 546             write (iout,'(a,i2,a,f10.5,1x,a)')
 547      &    ('b(',j,')=',b(j,i),aster(mask_fourier(j,i)),j=1,13)
 548 #endif
 549         enddo
 550 c Write a file with parameters for UNRES
 551         if (nparmset.eq.1) then
 552         open(88,file="fourier_opt.parm."//prefix(:ilen(prefix)),
 553      &     status="unknown")
 554         else
 555         open(88,file="fourier_opt.parm."//prefix(:ilen(prefix))//
 556      &     "_par"//liczba,
 557      &     status="unknown")
 558         endif
 559         write (88,'(i5,5x,a)') nloctyp,
 560      &                         "# Number of local interaction types"
 561 #ifdef NEWCORR
 562         write (88,'(40i2)') (itype2loc(i),i=1,ntyp)
 563         write (88,'(20i4)') (iloctyp(i),i=0,nloctyp)
 564         do i=0,nloctyp-1
 565           write (88,'(a)') restyp(iloctyp(i))
 566           do ii=1,3
 567             do j=1,2
 568               write (88,'(f20.10,4h  b1,i1,1h(,i1,1h))')
 569      &         bnew1(ii,j,i),j,ii
 570             enddo
 571           enddo
 572           do ii=1,3
 573             do j=1,2
 574               write (88,'(f20.10,4h  b2,i1,1h(,i1,1h))')
 575      &         bnew2(ii,j,i),j,ii
 576             enddo
 577           enddo
 578           do ii=1,3
 579             do j=1,2
 580               write (88,'(f20.10,4h  c1,i1,1h(,i1,1h))')
 581      &         2*ccnew(ii,j,i),j,ii
 582             enddo
 583           enddo
 584           do ii=1,3
 585             do j=1,2
 586               write (88,'(f20.10,4h  d1,i1,1h(,i1,1h))')
 587      &         2*ddnew(ii,j,i),j,ii
 588             enddo
 589           enddo
 590           do ii=1,2
 591             do j=1,2
 592               do k=1,2
 593               write (88,'(f20.10,3h  e,2i1,1h(,i1,1h))')
 594      &         eenew(ii,j,k,i),j,k,ii
 595               enddo
 596             enddo
 597           enddo
 598           do ii=1,3
 599               write (88,'(f20.10,5h  e0(,i1,1h))')
 600      &         e0new(ii,i),ii
 601           enddo
 602         enddo
 603 #else
 604         do i=1,nloctyp
 605           write (88,'(a)') restyp(iloctyp(i))
 606           do j=1,13
 607             write (88,'(f20.15)') b(j,i)
 608           enddo
 609         enddo
 610 #endif
 611         close(88)
 612       endif
 613       if (mod_elec) then
 614                 write (iout,'(/a)')
 615      &  'Electrostatic interaction constants (asterisk: optimized):'
 616         write (iout,'(1x,a,1x,a,7x,a,15x,a,15x,a,13x,a)')
 617      &            'IT','JT','EPP','RPP','ELPP6','ELPP3'
 618         do i=1,2
 619           do j=1,2
 620             write(iout,'(2i3,4(1pe15.4,1x,a1,1x))')i,j,
 621      &       epp(i,j),aster(mask_elec(i,j,1)),
 622      &       rpp(i,j),aster(mask_elec(i,j,2)),
 623      &       elpp6(i,j),aster(mask_elec(i,j,3)),
 624      &       elpp3(i,j),aster(mask_elec(i,j,4))
 625           enddo
 626         enddo
 627         write (iout,*)
 628         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
 629      &            'IT','JT','APP','BPP','AEL6','AEL3'
 630         do i=1,2
 631           do j=1,2
 632             write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
 633      &                    ael6(i,j),ael3(i,j)
 634           enddo
 635         enddo
 636 C Write electrostatic interaction parameters to a file
 637         if (nparmset.eq.1) then
 638         open(88,file="electr_opt.parm."//prefix(:ilen(prefix)),
 639      &   status="unknown")
 640         else
 641         open(88,file="electr_opt.parm."//prefix(:ilen(prefix))//
 642      &   "_par"//liczba,
 643      &   status="unknown")
 644         endif
 645         write(88,'(4f15.10,5x,a)')((epp(i,j),j=1,2),i=1,2),"! EPP"
 646         write(88,'(4f15.10,5x,a)')((rpp(i,j),j=1,2),i=1,2),"! RPP"
 647         write(88,'(4f15.10,5x,a)')((elpp6(i,j),j=1,2),i=1,2),"! ELPP6"
 648         write(88,'(4f15.10,5x,a)')((elpp3(i,j),j=1,2),i=1,2),"! ELPP3"
 649         close(88)
 650       endif
 651       if (mod_scp) then
 652         write (iout,*)
 653         write (iout,*)
 654      &  "Parameters of SC-p interactions (star: optimized):"
 655         write (iout,'(t14,a,t34,a)') "TYPE 1","TYPE 2"
 656         write (iout,'(8a)') "SC TYPE"," EPS_SCP  ","   R_SCP  ",
 657      &    " EPS_SCP  ","   R_SCP  "
 658         do i=1,20
 659           ast11 = aster(mask_scp(i,1,1))
 660           ast12 = aster(mask_scp(i,1,2))
 661           ast21 = aster(mask_scp(i,2,1))
 662           ast22 = aster(mask_scp(i,2,2))
 663           if (mask_scp(i,1,1).eq.0) ast11 = aster(mask_scp(0,1,1))
 664           if (mask_scp(i,1,2).eq.0) ast12 = aster(mask_scp(0,1,2))
 665           if (mask_scp(i,2,1).eq.0) ast21 = aster(mask_scp(0,2,1))
 666           if (mask_scp(i,2,2).eq.0) ast22 = aster(mask_scp(0,2,2))
 667           write (iout,'(2x,A3,2x,4(f8.3,1x,a1))') restyp(i),
 668      &     eps_scp(i,1),ast11,rscp(i,1),ast12,eps_scp(i,2),ast21,
 669      &     rscp(i,2),ast22
 670         enddo
 671 C Write SCp parameters to a file
 672         if (nparmset.eq.1) then
 673         open(88,file="scp_opt.parm."//prefix(:ilen(prefix)),
 674      &     status="unknown")
 675         else
 676         open(88,file="scp_opt.parm."//prefix(:ilen(prefix))//
 677      &     "_par"//liczba,
 678      &     status="unknown")
 679         endif
 680         do i=1,ntyp
 681           write(88,'(4f15.10,5x,a)') eps_scp(i,1),rscp(i,1),
 682      &      eps_scp(i,2),rscp(i,2),restyp(i)
 683         enddo
 684         close(88)
 685       endif
 686 50    format(bz,15(f8.5,1x,a1))
 687 60    format(bz,100f10.5)
 688       end
 689 c------------------------------------------------------------------------------
 690       subroutine write_zscore
 691       implicit none
 692       include "DIMENSIONS"
 693       include "DIMENSIONS.ZSCOPT"
 694       include "COMMON.ENERGIES"
 695       include "COMMON.CLASSES"
 696       include "COMMON.PROTNAME"
 697       include "COMMON.VMCPAR"
 698       include "COMMON.IOUNITS"
 699       include "COMMON.WEIGHTS"
 700       include "COMMON.WEIGHTDER"
 701       include "COMMON.COMPAR"
 702       include "COMMON.OPTIM"
 703       include "COMMON.THERMAL"
 704       integer i,j,k,iprot,ib
 705       character*6 namnat(5) /"angrms","Q","rmsd","rgy","sign"/
 706       integer ilen
 707       external ilen
 708       do iprot=1,nprot
 709
 710         write (iout,*)
 711         write (iout,'(a,2x,a)') "Protein",
 712      &    protname(iprot)(:ilen(protname(iprot)))
 713
 714         write (iout,*)
 715         write (iout,'(a)')"Maximum likelihood."
 716         write (iout,*)
 717         do ib=1,nbeta(1,iprot)
 718           write (iout,'(f8.1,f10.5)')
 719      &      1.0d0/(Rgas*betaT(ib,1,iprot)),sumlik(ib,iprot)
 720         enddo
 721         write (iout,*)
 722
 723         write (iout,*)
 724         write (iout,'(a)')"Specific heat."
 725         write (iout,*)
 726         write (iout,'(a,f10.5)') "Baseline",heatbase(iprot)
 727         do ib=1,nbeta(2,iprot)
 728           write (iout,'(f8.1,3f10.5)')
 729      &      1.0d0/(Rgas*betaT(ib,2,iprot)),
 730      &      zvtot(ib,iprot),target_cv(ib,iprot),weiCv(ib,iprot)
 731         enddo
 732
 733         if (natlike(iprot).gt.0) write (iout,'(a)')"Nativelikeness."
 734         do i=1,natlike(iprot)
 735           write (iout,*) "Property",i,numnat(i,iprot)
 736           if (natdim(i,iprot).eq.1)  then
 737             do ib=1,nbeta(i+2,iprot)
 738               write (iout,'(f7.2,2f10.5)')
 739      &         1.0d0/(Rgas*betaT(ib,i+2,iprot)),nuave(1,ib,i,iprot),
 740      &         nuexp(1,ib,i,iprot),weinat(1,ib,i,iprot)
 741             enddo
 742           else
 743             do ib=1,nbeta(i+2,iprot)
 744               write (iout,'("Temperature",f7.2," NATDIM",i5)')
 745      &        1.0d0/(Rgas*betaT(ib,i+2,iprot)),natdim(i,iprot)
 746               do k=1,natdim(i,iprot)
 747               write (iout,'(3f10.5)') nuave(k,ib,i,iprot),
 748      &         nuexp(k,ib,i,iprot),weinat(k,ib,i,iprot)
 749               enddo
 750             enddo
 751           endif
 752         enddo
 753
 754       enddo
 755       write (iout,*)
 756       return
 757       end
 758 c------------------------------------------------------------------------------
 759       subroutine write_prot(ii,przed)
 760       implicit none
 761       include "DIMENSIONS.ZSCOPT"
 762       include "COMMON.ENERGIES"
 763       include "COMMON.CLASSES"
 764       include "COMMON.PROTNAME"
 765       include "COMMON.VMCPAR"
 766       include "COMMON.IOUNITS"
 767       include "COMMON.COMPAR"
 768       integer i,ii,j
 769       integer ilen
 770       external ilen
 771       character*(*) przed
 772       character*64 nazwa
 773       nazwa=przed(:ilen(przed))//'.'//prefix(:ilen(prefix))
 774      &   //'.'//protname(ii)
 775       open(unit=istat,file=nazwa)
 776       do i=1,ntot_work(ii)
 777         write(istat,'(i10,3f15.3,20e15.5)')
 778      &    i,e_total(i,ii),eini(i,ii),entfac(i,ii),
 779      &    (Ptab(i,j,ii),j=1,nbeta(1,ii))
 780       enddo
 781       close(istat)
 782       return
 783       end