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