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