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