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