source/maxlik/src_FPy/maxlikopt.F

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