source/wham/src/enecalc1.F

   1       subroutine enecalc(islice,*)
   2       implicit none
   3       include "DIMENSIONS"
   4       include "DIMENSIONS.ZSCOPT"
   5       include "DIMENSIONS.FREE"
   6 #ifdef MPI
   7       include "mpif.h"
   8       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
   9       include "COMMON.MPI"
  10 #endif
  11       include "COMMON.CHAIN"
  12       include "COMMON.IOUNITS"
  13       include "COMMON.PROTFILES"
  14       include "COMMON.NAMES"
  15       include "COMMON.VAR"
  16       include "COMMON.SBRIDGE"
  17       include "COMMON.GEO"
  18       include "COMMON.FFIELD"
  19       include "COMMON.ENEPS"
  20       include "COMMON.LOCAL"
  21       include "COMMON.WEIGHTS"
  22       include "COMMON.INTERACT"
  23       include "COMMON.FREE"
  24       include "COMMON.ENERGIES"
  25       include "COMMON.CONTROL"
  26       include "COMMON.TORCNSTR"
  27       character*64 nazwa
  28       character*80 bxname
  29       character*3 liczba
  30       double precision qwolynes
  31       external qwolynes
  32       integer errmsg_count,maxerrmsg_count /100/
  33       double precision rmsnat,gyrate
  34       external rmsnat,gyrate
  35       double precision tole /1.0d-1/
  36       integer i,itj,ii,iii,j,k,l,licz,scme,itmp
  37       integer ires
  38       integer ir,ib,ipar,iparm
  39       integer iscor,islice
  40       real*4 csingle(3,maxres2)
  41       double precision energ
  42       integer ilen,iroof
  43       external ilen,iroof
  44       double precision energia(0:max_ene),rmsdev,efree,eini
  45       double precision fT(6),quot,quotl,kfacl,kfac /2.4d0/,T0 /3.0d2/
  46       double precision tt
  47       integer snk_p(MaxR,MaxT_h,Max_parm)
  48       logical lerr
  49       character*64 bprotfile_temp
  50       integer scount_t(0:maxprocs-1)
  51       call opentmp(islice,ientout,bprotfile_temp)
  52       iii=0
  53       ii=0
  54       errmsg_count=0
  55       write (iout,*) "enecalc: nparmset ",nparmset
  56 #ifdef MPI
  57       do iparm=1,nParmSet
  58         do ib=1,nT_h(iparm)
  59           do i=1,nR(ib,iparm)
  60             snk_p(i,ib,iparm)=0
  61           enddo
  62         enddo
  63       enddo
  64       do i=indstart(me1),indend(me1)
  65 #else
  66       do iparm=1,nParmSet
  67         do ib=1,nT_h(iparm)
  68           do i=1,nR(ib,iparm)
  69             snk(i,ib,iparm)=0
  70           enddo
  71         enddo
  72       enddo
  73       do i=1,ntot
  74 #endif
  75         read(ientout,rec=i,err=101)
  76      &    ((csingle(l,k),l=1,3),k=1,nres),
  77      &    ((csingle(l,k+nres),l=1,3),k=nnt,nct),
  78      &    nss,(ihpb(k),jhpb(k),k=1,nss),
  79      &    eini,efree,rmsdev,(q(j,iii+1),j=1,nQ),iR,ib,ipar
  80 cc       write(iout,*), 'NAWEJ',i,eini
  81          if (indpdb.gt.0) then
  82            do k=1,nres
  83              do l=1,3
  84                c(l,k)=csingle(l,k)
  85              enddo
  86            enddo
  87            do k=nnt,nct
  88              do l=1,3
  89                c(l,k+nres)=csingle(l,k+nres)
  90              enddo
  91            enddo
  92            q(nQ+1,iii+1)=rmsnat(iii+1)
  93          endif
  94          q(nQ+2,iii+1)=gyrate(iii+1)
  95 c        fT=T0*beta_h(ib,ipar)*1.987D-3
  96 c        ft=2.0d0/(1.0d0+1.0d0/(T0*beta_h(ib,ipar)*1.987D-3))
  97         if (rescale_mode.eq.1) then
  98           quot=1.0d0/(T0*beta_h(ib,ipar)*1.987D-3)
  99 #if defined(FUNCTH)
 100           tt = 1.0d0/(beta_h(ib,ipar)*1.987D-3)
 101           ft(6)=(320.0+80.0*dtanh((tt-320.0)/80.0))/320.0
 102 #elif defined(FUNCT)
 103           ft(6)=quot
 104 #else
 105           ft(6)=1.0d0
 106 #endif
 107           quotl=1.0d0
 108           kfacl=1.0d0
 109           do l=1,5
 110             quotl=quotl*quot
 111             kfacl=kfacl*kfac
 112             fT(l)=kfacl/(kfacl-1.0d0+quotl)
 113           enddo
 114         else if (rescale_mode.eq.2) then
 115           quot=1.0d0/(T0*beta_h(ib,ipar)*1.987D-3)
 116 #if defined(FUNCTH)
 117           tt = 1.0d0/(beta_h(ib,ipar)*1.987D-3)
 118           ft(6)=(320.0+80.0*dtanh((tt-320.0)/80.0))/320.0
 119 #elif defined(FUNCT)
 120           ft(6)=quot
 121 #else
 122           ft(6)=1.0d0
 123 #endif
 124           quotl=1.0d0
 125           do l=1,5
 126             quotl=quotl*quot
 127             fT(l)=1.12692801104297249644d0/
 128      &         dlog(dexp(quotl)+dexp(-quotl))
 129           enddo
 130         else if (rescale_mode.eq.0) then
 131           do l=1,5
 132             fT(l)=1.0d0
 133           enddo
 134         else
 135           write (iout,*) "Error in ECECALC: wrong RESCALE_MODE",
 136      &     rescale_mode
 137           call flush(iout)
 138           return1
 139         endif
 140
 141 c        write (iout,*) "T",1.0d0/(beta_h(ib,ipar)*1.987D-3)," T0",T0,
 142 c     &   " kfac",kfac,"quot",quot," fT",fT
 143         do j=1,2*nres
 144           do k=1,3
 145             c(k,j)=csingle(k,j)
 146           enddo
 147         enddo
 148         call int_from_cart1(.false.)
 149         ii=ii+1
 150         do iparm=1,nparmset
 151
 152         call restore_parm(iparm)
 153 #ifdef DEBUG
 154             write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
 155      &      wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
 156      &      wtor_d,wsccor,wbond
 157 #endif
 158         call etotal(energia(0),fT)
 159 #ifdef DEBUG
 160         write (iout,*) "Conformation",i
 161         call enerprint(energia(0),fT)
 162 c        write (iout,'(2i5,21f8.2)') i,iparm,(energia(k),k=1,21)
 163 c        write (iout,*) "ftors",ftors
 164       write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
 165      & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
 166       do ires=1,nres
 167         write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
 168      &    restyp(itype(ires)),ires,(c(j,ires),j=1,3),
 169      &    (c(j,ires+nres),j=1,3)
 170       enddo
 171         call intout
 172 #endif
 173         if (energia(0).ge.1.0d20) then
 174           write (iout,*) "NaNs detected in some of the energy",
 175      &     " components for conformation",ii+1
 176           write (iout,*) "The Cartesian geometry is:"
 177           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 178           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 179           write (iout,*) "The internal geometry is:"
 180 c          call intout
 181 c        call pdbout(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
 182           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 183           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 184           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 185           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 186           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 187           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 188           write (iout,*) "The components of the energy are:"
 189           call enerprint(energia(0),fT)
 190           write (iout,*)
 191      &      "This conformation WILL NOT be added to the database."
 192           call flush(iout)
 193           goto 121
 194         else
 195 #ifdef DEBUG
 196           if (ipar.eq.iparm) write (iout,*) i,iparm,
 197      &      1.0d0/(beta_h(ib,ipar)*1.987D-3),eini,energia(0)
 198 #endif
 199           if (ipar.eq.iparm .and. einicheck.gt.0 .and.
 200      &      dabs(eini-energia(0)).gt.tole) then
 201             if (errmsg_count.le.maxerrmsg_count) then
 202               write (iout,'(2a,2e15.5,a,2i8,a,f8.1)')
 203      &         "Warning: energy differs remarkably from ",
 204      &         " the value read in: ",energia(0),eini," point",
 205      &         iii+1,indstart(me1)+iii," T",
 206      &         1.0d0/(1.987D-3*beta_h(ib,ipar))
 207 #ifdef DEBUG
 208              call enerprint(energia(0),fT)
 209              write (iout,'(4f10.5,2i5)') 0.0,energia(0),0.0,
 210      &       1.0d0/(beta_h(ib,ipar)*1.987D-3),
 211      &       0,0
 212              write(iout,'(8f10.5)')
 213      &       ((c(l,k),l=1,3),k=1,nres),
 214      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
 215              itmp=ipdb
 216              ipdb=iout
 217              call pdbout(iii+1,beta_h(ib,ipar),
 218      &                   eini,energia(0),0.0d0,rmsdev)
 219              write (iout,*)
 220              ipdb=itmp
 221 #endif
 222               errmsg_count=errmsg_count+1
 223               if (errmsg_count.gt.maxerrmsg_count)
 224      &          write (iout,*) "Too many warning messages"
 225               if (einicheck.gt.1) then
 226                 write (iout,*) "Calculation stopped."
 227                 call flush(iout)
 228 #ifdef MPI
 229                 call MPI_Abort(WHAM_COMM,IERROR,ERRCODE)
 230 #endif
 231                 call flush(iout)
 232                 return1
 233               endif
 234             endif
 235           endif
 236           potE(iii+1,iparm)=energia(0)
 237           do k=1,max_ene
 238             enetb(k,iii+1,iparm)=energia(k)
 239           enddo
 240 c          write(iout,*) "iCHUJ TU STRZELI",i,iii,entfac(i)
 241 c          call enerprint(energia(0),fT)
 242 #ifdef DEBUG
 243           write (iout,'(2i5,f10.1,3e15.5)') i,iii,
 244      &     1.0d0/(beta_h(ib,ipar)*1.987D-3),energia(0),eini,efree
 245           call enerprint(energia(0),fT)
 246           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 247           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 248           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 249           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 250           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 251           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 252           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 253           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 254           write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
 255           write (iout,'(8f10.5)') (q(k,iii+1),k=1,nQ)
 256           write (iout,'(f10.5,i10)') rmsdev,iscor
 257           call enerprint(energia(0),fT)
 258         write(liczba,'(bz,i3.3)') me
 259         nazwa="test"//liczba//".pdb"
 260         write (iout,*) "pdb file",nazwa
 261         open (ipdb,file=nazwa,position="append")
 262         call pdbout(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
 263         close(ipdb)
 264 #endif
 265         endif
 266
 267         enddo ! iparm
 268
 269         iii=iii+1
 270         if (q(1,iii).le.0.0d0 .and. indpdb.gt.0) q(1,iii)=qwolynes(0,0)
 271         write (ientout,rec=iii)
 272      &   ((csingle(l,k),l=1,3),k=1,nres),
 273      &   ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 274      &   nss,(ihpb(k),jhpb(k),k=1,nss),
 275      &   potE(iii,ipar),efree,rmsdev,(q(k,iii),k=1,nQ),iR,ib,ipar
 276 c        write (iout,'(2i5,2e15.5)') ii,iii,potE(iii,ipar),efree
 277 #ifdef MPI
 278         if (separate_parset) then
 279           snk_p(iR,ib,1)=snk_p(iR,ib,1)+1
 280         else
 281           snk_p(iR,ib,ipar)=snk_p(iR,ib,ipar)+1
 282         endif
 283 c        write (iout,*) "iii",iii," iR",iR," ib",ib," ipar",ipar,
 284 c     &   " snk",snk_p(iR,ib,ipar)
 285 #else
 286         snk(iR,ib,ipar,islice)=snk(iR,ib,ipar,islice)+1
 287 #endif
 288   121   continue
 289       enddo
 290 #ifdef MPI
 291       scount(me)=iii
 292       write (iout,*) "Me",me," scount",scount(me)
 293       call flush(iout)
 294 c  Master gathers updated numbers of conformations written by all procs.
 295       scme = scount(me)
 296       call MPI_AllGather( scount(me), 1, MPI_INTEGER, scount_t(0), 1,
 297      &  MPI_INTEGER, WHAM_COMM, IERROR)
 298       do k=0,nprocs-1
 299         scount(k) = scount_t(k)
 300       enddo
 301       indstart(0)=1
 302       indend(0)=scount(0)
 303       do i=1, Nprocs-1
 304         indstart(i)=indend(i-1)+1
 305         indend(i)=indstart(i)+scount(i)-1
 306       enddo
 307       write (iout,*)
 308       write (iout,*) "Revised conformation counts"
 309       do i=0,nprocs1-1
 310         write (iout,'(a,i5,a,i7,a,i7,a,i7)')
 311      &    "Processor",i," indstart",indstart(i),
 312      &    " indend",indend(i)," count",scount(i)
 313       enddo
 314       call flush(iout)
 315       call MPI_AllReduce(snk_p(1,1,1),snk(1,1,1,islice),
 316      &  MaxR*MaxT_h*nParmSet,
 317      &  MPI_INTEGER,MPI_SUM,WHAM_COMM,IERROR)
 318 #endif
 319       stot(islice)=0
 320       do iparm=1,nParmSet
 321         do ib=1,nT_h(iparm)
 322           do i=1,nR(ib,iparm)
 323             stot(islice)=stot(islice)+snk(i,ib,iparm,islice)
 324           enddo
 325         enddo
 326       enddo
 327       write (iout,*) "Revised SNK"
 328       do iparm=1,nParmSet
 329         do ib=1,nT_h(iparm)
 330           write (iout,'("Param",i3," Temp",f6.1,3x,32i8)')
 331      &     iparm,1.0d0/(1.987D-3*beta_h(ib,iparm)),
 332      &     (snk(i,ib,iparm,islice),i=1,nR(ib,iparm))
 333           write (iout,*) "snk_p",(snk_p(i,ib,iparm),i=1,nR(ib,iparm))
 334         enddo
 335       enddo
 336       write (iout,'("Total",i10)') stot(islice)
 337       call flush(iout)
 338       return
 339   101 write (iout,*) "Error in scratchfile."
 340       call flush(iout)
 341       return1
 342       end
 343 c------------------------------------------------------------------------------
 344       subroutine write_dbase(islice,*)
 345       implicit none
 346       include "DIMENSIONS"
 347       include "DIMENSIONS.ZSCOPT"
 348       include "DIMENSIONS.FREE"
 349       include "DIMENSIONS.COMPAR"
 350 #ifdef MPI
 351       include "mpif.h"
 352       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 353       include "COMMON.MPI"
 354 #endif
 355       include "COMMON.CONTROL"
 356       include "COMMON.CHAIN"
 357       include "COMMON.IOUNITS"
 358       include "COMMON.PROTFILES"
 359       include "COMMON.NAMES"
 360       include "COMMON.VAR"
 361       include "COMMON.SBRIDGE"
 362       include "COMMON.GEO"
 363       include "COMMON.FFIELD"
 364       include "COMMON.ENEPS"
 365       include "COMMON.LOCAL"
 366       include "COMMON.WEIGHTS"
 367       include "COMMON.INTERACT"
 368       include "COMMON.FREE"
 369       include "COMMON.ENERGIES"
 370       include "COMMON.COMPAR"
 371       include "COMMON.PROT"
 372       character*64 nazwa
 373       character*80 bxname,cxname
 374       character*64 bprotfile_temp
 375       character*3 liczba,licz
 376       character*2 licz2
 377       integer i,itj,ii,iii,j,k,l
 378       integer ixdrf,iret
 379       integer iscor,islice
 380       double precision rmsdev,efree,eini
 381       real*4 csingle(3,maxres2)
 382       double precision energ
 383       integer ilen,iroof
 384       external ilen,iroof
 385       integer ir,ib,iparm
 386       write (licz2,'(bz,i2.2)') islice
 387       call opentmp(islice,ientout,bprotfile_temp)
 388       write (iout,*) "bprotfile_temp ",bprotfile_temp
 389       call flush(iout)
 390       if (.not.bxfile .and. .not. cxfile .and. indpdb.eq.0
 391      &   .and. ensembles.eq.0) then
 392         close(ientout,status="delete")
 393         return
 394       endif
 395 #ifdef MPI
 396       write (liczba,'(bz,i3.3)') me
 397       if (bxfile .or. cxfile .or. ensembles.gt.0) then
 398         if (.not.separate_parset) then
 399           bxname = prefix(:ilen(prefix))//liczba//".bx"
 400         else
 401           write (licz,'(bz,i3.3)') myparm
 402           bxname = prefix(:ilen(prefix))//liczba//"_par"//licz//".bx"
 403         endif
 404         open (ientin,file=bxname,status="unknown",
 405      &    form="unformatted",access="direct",recl=lenrec1)
 406       endif
 407 #else
 408       if (bxfile .or. cxfile .or. ensembles.gt.0) then
 409         if (nslice.eq.1) then
 410           bxname = prefix(:ilen(prefix))//".bx"
 411         else
 412           bxname = prefix(:ilen(prefix))//
 413      &           "_slice_"//licz2//".bx"
 414         endif
 415         open (ientin,file=bxname,status="unknown",
 416      &    form="unformatted",access="direct",recl=lenrec1)
 417         write (iout,*) "Calculating energies; writing geometry",
 418      & " and energy components to ",bxname(:ilen(bxname))
 419       endif
 420 #if (defined(AIX) && !defined(JUBL))
 421         call xdrfopen_(ixdrf,cxname, "w", iret)
 422 #else
 423         call xdrfopen(ixdrf,cxname, "w", iret)
 424 #endif
 425         if (iret.eq.0) then
 426           write (iout,*) "Error opening cxfile ",cxname(:ilen(cxname))
 427           cxfile=.fale.
 428         endif
 429       endif
 430 #endif
 431       if (indpdb.gt.0) then
 432         if (nslice.eq.1) then
 433 #ifdef MPI
 434          if (.not.separate_parset) then
 435            statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))
 436      &       //liczba//'.stat'
 437          else
 438            write (licz,'(bz,i3.3)') myparm
 439            statname=prefix(:ilen(prefix))//'_par'//licz//'_'//
 440      &      pot(:ilen(pot))//liczba//'.stat'
 441          endif
 442
 443 #else
 444           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))//'.stat'
 445 #endif
 446         else
 447 #ifdef MPI
 448          if (.not.separate_parset) then
 449           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))//
 450      &      "_slice_"//licz2//liczba//'.stat'
 451          else
 452           write (licz,'(bz,i3.3)') myparm
 453           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))//
 454      &      '_par'//licz//"_slice_"//licz2//liczba//'.stat'
 455          endif
 456 #else
 457           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))
 458      &      //"_slice_"//licz2//'.stat'
 459 #endif
 460         endif
 461         open(istat,file=statname,status="unknown")
 462       endif
 463
 464 #ifdef MPI
 465       do i=1,scount(me)
 466 #else
 467       do i=1,ntot(islice)
 468 #endif
 469 cc        if (dyn_ss) then
 470 cc        read(ientout,rec=i,err=101)
 471 cc     &    ((csingle(l,k),l=1,3),k=1,nres),
 472 cc     &    ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 473 cc     &    nss,(idssb(k),jdssb(k),k=1,nss),
 474 cc     &    eini,efree,rmsdev,(q(k,i),k=1,nQ),iR,ib,iparm
 475 cc        idssb(k)=idssb(k)-nres
 476 cc        jdssb(k)=jdssb(k)-nres
 477 cc        else
 478         read(ientout,rec=i,err=101)
 479      &    ((csingle(l,k),l=1,3),k=1,nres),
 480      &    ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 481      &    nss,(ihpb(k),jhpb(k),k=1,nss),
 482      &    eini,efree,rmsdev,(q(k,i),k=1,nQ),iR,ib,iparm
 483 cc         endif
 484 cc        write (iout,*) 'CC', iR,ib,iparm,eini,efree
 485         do j=1,2*nres
 486           do k=1,3
 487             c(k,j)=csingle(k,j)
 488           enddo
 489         enddo
 490         call int_from_cart1(.false.)
 491         iscore=0
 492         if (indpdb.gt.0) then
 493           call conf_compar(i,.false.,.true.)
 494         endif
 495 c        if (dyn_ss) then
 496         if (bxfile .or.cxfile .or. ensembles.gt.0) write (ientin,rec=i)
 497      &    ((csingle(l,k),l=1,3),k=1,nres),
 498      &    ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 499      &    nss,(ihpb(k),jhpb(k),k=1,nss),
 500 c     &    potE(i,iparm),-entfac(i),rms_nat,iscore
 501      &    potE(i,nparmset),-entfac(i),rms_nat,iscore
 502 c        else
 503           if (bxfile .or.cxfile .or. ensembles.gt.0) write
 504      &     (ientin,rec=i)
 505      &    ((csingle(l,k),l=1,3),k=1,nres),
 506      &    ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 507      &    nss,(ihpb(k),jhpb(k),k=1,nss),
 508 c     &    potE(i,iparm),-entfac(i),rms_nat,iscore
 509      &    potE(i,nparmset),-entfac(i),rms_nat,iscore
 510 c        endif
 511         write (iout,'(2i5,3e15.5)') i,me,potE(i,iparm),-entfac(i)
 512 #ifndef MPI
 513         if (cxfile) call cxwrite(ixdrf,csingle,potE(i,nparmset),
 514      &    -entfac(i),rms_nat,iscore)
 515 #endif
 516       enddo
 517       close(ientout,status="delete")
 518       close(istat)
 519       if (bxfile .or. cxfile .or. ensembles.gt.0) close(ientin)
 520 #ifdef MPI
 521       call MPI_Barrier(WHAM_COMM,IERROR)
 522       if (me.ne.Master .or. .not.bxfile .and. .not. cxfile
 523      &   .and. ensembles.eq.0) return
 524       write (iout,*)
 525       if (bxfile .or. ensembles.gt.0) then
 526         if (nslice.eq.1) then
 527           if (.not.separate_parset) then
 528             bxname = prefix(:ilen(prefix))//".bx"
 529           else
 530             write (licz,'(bz,i3.3)') myparm
 531             bxname = prefix(:ilen(prefix))//"_par"//licz//".bx"
 532           endif
 533         else
 534           if (.not.separate_parset) then
 535             bxname = prefix(:ilen(prefix))//"_slice_"//licz2//".bx"
 536           else
 537             write (licz,'(bz,i3.3)') myparm
 538             bxname = prefix(:ilen(prefix))//"par_"//licz//
 539      &        "_slice_"//licz2//".bx"
 540           endif
 541         endif
 542         open (ientout,file=bxname,status="unknown",
 543      &      form="unformatted",access="direct",recl=lenrec1)
 544         write (iout,*) "Master is creating binary database ",
 545      &   bxname(:ilen(bxname))
 546       endif
 547       if (cxfile) then
 548         if (nslice.eq.1) then
 549           if (.not.separate_parset) then
 550             cxname = prefix(:ilen(prefix))//".cx"
 551           else
 552             cxname = prefix(:ilen(prefix))//"_par"//licz//".cx"
 553           endif
 554         else
 555           if (.not.separate_parset) then
 556             cxname = prefix(:ilen(prefix))//
 557      &             "_slice_"//licz2//".cx"
 558           else
 559             cxname = prefix(:ilen(prefix))//"_par"//licz//
 560      &             "_slice_"//licz2//".cx"
 561           endif
 562         endif
 563 #if (defined(AIX) && !defined(JUBL))
 564         call xdrfopen_(ixdrf,cxname, "w", iret)
 565 #else
 566         call xdrfopen(ixdrf,cxname, "w", iret)
 567 #endif
 568         if (iret.eq.0) then
 569           write (iout,*) "Error opening cxfile ",cxname(:ilen(cxname))
 570           cxfile=.false.
 571         endif
 572       endif
 573       do j=0,nprocs-1
 574         write (liczba,'(bz,i3.3)') j
 575         if (separate_parset) then
 576           write (licz,'(bz,i3.3)') myparm
 577           bxname = prefix(:ilen(prefix))//liczba//"_par"//licz//".bx"
 578         else
 579           bxname = prefix(:ilen(prefix))//liczba//".bx"
 580         endif
 581         open (ientin,file=bxname,status="unknown",
 582      &    form="unformatted",access="direct",recl=lenrec1)
 583         write (iout,*) "Master is reading conformations from ",
 584      &   bxname(:ilen(bxname))
 585         iii = 0
 586 c        write (iout,*) "j",j," indstart",indstart(j)," indend",indend(j)
 587 c        call flush(iout)
 588         do i=indstart(j),indend(j)
 589           iii = iii+1
 590 cc          if (dyn_ss) then
 591 cc          read(ientin,rec=iii,err=101)
 592 cc     &      ((csingle(l,k),l=1,3),k=1,nres),
 593 cc     &      ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 594 cc     &      nss,(idssb(k),jdssb(k),k=1,nss),
 595 cc    &      eini,efree,rmsdev,iscor
 596 cc          idssb(k)=idssb(k)-nres
 597 cc          jdssb(k)=jdssb(k)-nres
 598 cc          else
 599           read(ientin,rec=iii,err=101)
 600      &      ((csingle(l,k),l=1,3),k=1,nres),
 601      &      ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 602      &      nss,(ihpb(k),jhpb(k),k=1,nss),
 603      &      eini,efree,rmsdev,iscor
 604 cc          endif
 605           if (bxfile .or. ensembles.gt.0) then
 606 cc          if (dyn_ss) then
 607 cc            write (ientout,rec=i)
 608 cc     &        ((csingle(l,k),l=1,3),k=1,nres),
 609 cc     &        ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 610 cc     &        nss,(idssb(k)+nres,jdssb(k)+nres,k=1,nss),
 611 cc    &        eini,efree,rmsdev,iscor
 612 cc           else
 613                         write (ientout,rec=i)
 614      &        ((csingle(l,k),l=1,3),k=1,nres),
 615      &        ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 616      &        nss,(ihpb(k),jhpb(k),k=1,nss),
 617      &        eini,efree,rmsdev,iscor
 618 cc           write(iout,*) "W poszukiwaniu zlotych galotow"
 619 cc           write(iout,*) "efree=",efree,iii
 620 cc           endif
 621           endif
 622           if(cxfile)call cxwrite(ixdrf,csingle,eini,efree,rmsdev,iscor)
 623 #ifdef DEBUG
 624           do k=1,2*nres
 625             do l=1,3
 626               c(l,k)=csingle(l,k)
 627             enddo
 628           enddo
 629           call int_from_cart1(.false.)
 630           write (iout,'(2i5,3e15.5)') i,iii,eini,efree
 631           write (iout,*) "The Cartesian geometry is:"
 632           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 633           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 634           write (iout,*) "The internal geometry is:"
 635           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 636           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 637           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 638           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 639           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 640           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 641           write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
 642           write (iout,'(f10.5,i5)') rmsdev,iscor
 643 #endif
 644         enddo ! i
 645         write (iout,*) iii," conformations (from",indstart(j)," to",
 646      &   indend(j),") read from ",
 647      &   bxname(:ilen(bxname))
 648         close (ientin,status="delete")
 649       enddo ! j
 650       if (bxfile .or. cxfile .or. ensembles.gt.0) close (ientout)
 651 #if (defined(AIX) && !defined(JUBL))
 652       if (cxfile) call xdrfclose_(ixdrf,cxname,iret)
 653 #else
 654       if (cxfile) call xdrfclose(ixdrf,cxname,iret)
 655 #endif
 656 #endif
 657       return
 658   101 write (iout,*) "Error in scratchfile."
 659       call flush(iout)
 660       return1
 661       end
 662 c-------------------------------------------------------------------------------
 663       subroutine cxwrite(ixdrf,csingle,eini,efree,rmsdev,iscor)
 664       implicit none
 665       include "DIMENSIONS"
 666       include "DIMENSIONS.ZSCOPT"
 667       include "DIMENSIONS.FREE"
 668       include "DIMENSIONS.COMPAR"
 669 #ifdef MPI
 670       include "mpif.h"
 671       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 672       include "COMMON.MPI"
 673 #endif
 674       include "COMMON.CONTROL"
 675       include "COMMON.CHAIN"
 676       include "COMMON.IOUNITS"
 677       include "COMMON.PROTFILES"
 678       include "COMMON.NAMES"
 679       include "COMMON.VAR"
 680       include "COMMON.SBRIDGE"
 681       include "COMMON.GEO"
 682       include "COMMON.FFIELD"
 683       include "COMMON.ENEPS"
 684       include "COMMON.LOCAL"
 685       include "COMMON.WEIGHTS"
 686       include "COMMON.INTERACT"
 687       include "COMMON.FREE"
 688       include "COMMON.ENERGIES"
 689       include "COMMON.COMPAR"
 690       include "COMMON.PROT"
 691       integer i,j,itmp,iscor,iret,ixdrf
 692       double precision rmsdev,efree,eini
 693       real*4 csingle(3,maxres2),xoord(3,maxres2+2)
 694       real*4 prec
 695
 696 c      write (iout,*) "cxwrite"
 697 c      call flush(iout)
 698       prec=10000.0
 699       do i=1,nres
 700        do j=1,3
 701         xoord(j,i)=csingle(j,i)
 702        enddo
 703       enddo
 704       do i=nnt,nct
 705        do j=1,3
 706         xoord(j,nres+i-nnt+1)=csingle(j,i+nres)
 707        enddo
 708       enddo
 709
 710       itmp=nres+nct-nnt+1
 711
 712 c      write (iout,*) "itmp",itmp
 713 c      call flush(iout)
 714 c       write (iout,*) "CNZ",eini,dyn_ss
 715 #if (defined(AIX) && !defined(JUBL))
 716       call xdrf3dfcoord_(ixdrf, xoord, itmp, prec, iret)
 717
 718 c      write (iout,*) "xdrf3dfcoord"
 719 c      call flush(iout)
 720       call xdrfint_(ixdrf, nss, iret)
 721       do j=1,nss
 722 cc       if (dyn_ss) then
 723 cc       call xdrfint_(ixdrf, idssb(j)+nres, iret)
 724 cc        call xdrfint_(ixdrf, jdssb(j)+nres, iret)
 725 cc       else
 726         call xdrfint_(ixdrf, ihpb(j), iret)
 727         call xdrfint_(ixdrf, jhpb(j), iret)
 728 cc       endif
 729       enddo
 730       call xdrffloat_(ixdrf,real(eini),iret)
 731       call xdrffloat_(ixdrf,real(efree),iret)
 732       call xdrffloat_(ixdrf,real(rmsdev),iret)
 733       call xdrfint_(ixdrf,iscor,iret)
 734 #else
 735       call xdrf3dfcoord(ixdrf, xoord, itmp, prec, iret)
 736
 737       call xdrfint(ixdrf, nss, iret)
 738       do j=1,nss
 739 cc       if (dyn_ss) then
 740 cc        call xdrfint(ixdrf, idssb(j), iret)
 741 cc        call xdrfint(ixdrf, jdssb(j), iret)
 742 cc        idssb(j)=idssb(j)-nres
 743 cc       jdssb(j)=jdssb(j)-nres
 744 cc       else
 745         call xdrfint(ixdrf, ihpb(j), iret)
 746         call xdrfint(ixdrf, jhpb(j), iret)
 747 cc       endif
 748       enddo
 749       call xdrffloat(ixdrf,real(eini),iret)
 750       call xdrffloat(ixdrf,real(efree),iret)
 751       call xdrffloat(ixdrf,real(rmsdev),iret)
 752       call xdrfint(ixdrf,iscor,iret)
 753 #endif
 754
 755       return
 756       end
 757 c------------------------------------------------------------------------------
 758       logical function conf_check(ii,iprint)
 759       implicit none
 760       include "DIMENSIONS"
 761       include "DIMENSIONS.ZSCOPT"
 762       include "DIMENSIONS.FREE"
 763 #ifdef MPI
 764       include "mpif.h"
 765       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 766       include "COMMON.MPI"
 767 #endif
 768       include "COMMON.CHAIN"
 769       include "COMMON.IOUNITS"
 770       include "COMMON.PROTFILES"
 771       include "COMMON.NAMES"
 772       include "COMMON.VAR"
 773       include "COMMON.SBRIDGE"
 774       include "COMMON.GEO"
 775       include "COMMON.FFIELD"
 776       include "COMMON.ENEPS"
 777       include "COMMON.LOCAL"
 778       include "COMMON.WEIGHTS"
 779       include "COMMON.INTERACT"
 780       include "COMMON.FREE"
 781       include "COMMON.ENERGIES"
 782       include "COMMON.CONTROL"
 783       include "COMMON.TORCNSTR"
 784       integer j,k,l,ii,itj,iprint
 785       if (.not.check_conf) then
 786         conf_check=.true.
 787         return
 788       endif
 789       call int_from_cart1(.false.)
 790       do j=nnt+1,nct
 791         if (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0) then
 792           if (iprint.gt.0)
 793      &    write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),
 794      &      " for conformation",ii
 795           if (iprint.gt.1) then
 796             write (iout,*) "The Cartesian geometry is:"
 797             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 798             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 799             write (iout,*) "The internal geometry is:"
 800             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 801             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 802             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 803             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 804             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 805             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 806           endif
 807           if (iprint.gt.0) write (iout,*)
 808      &      "This conformation WILL NOT be added to the database."
 809           conf_check=.false.
 810           return
 811         endif
 812       enddo
 813       do j=nnt,nct
 814         itj=itype(j)
 815         if (itype(j).ne.10 .and. (vbld(nres+j)-dsc(itj)).gt.2.0d0) then
 816           if (iprint.gt.0)
 817      &    write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j),
 818      &     " for conformation",ii
 819           if (iprint.gt.1) then
 820             write (iout,*) "The Cartesian geometry is:"
 821             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 822             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 823             write (iout,*) "The internal geometry is:"
 824             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 825             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 826             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 827             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 828             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 829             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 830           endif
 831           if (iprint.gt.0) write (iout,*)
 832      &      "This conformation WILL NOT be added to the database."
 833           conf_check=.false.
 834           return
 835         endif
 836       enddo
 837       do j=3,nres
 838         if (theta(j).le.0.0d0) then
 839           if (iprint.gt.0)
 840      &    write (iout,*) "Zero theta angle(s) in conformation",ii
 841           if (iprint.gt.1) then
 842             write (iout,*) "The Cartesian geometry is:"
 843             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 844             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 845             write (iout,*) "The internal geometry is:"
 846             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 847             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 848             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 849             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 850             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 851             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 852           endif
 853           if (iprint.gt.0) write (iout,*)
 854      &      "This conformation WILL NOT be added to the database."
 855           conf_check=.false.
 856           return
 857         endif
 858         if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
 859       enddo
 860       conf_check=.true.
 861 c      write (iout,*) "conf_check passed",ii
 862       return
 863       end