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