source/wham/src-M/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       integer scount_(0:MaxProcs)
  49       logical lerr
  50       character*64 bprotfile_temp
  51       call opentmp(islice,ientout,bprotfile_temp)
  52       iii=0
  53       ii=0
  54       errmsg_count=0
  55       write (iout,*) "enecalc: nparmset ",nparmset
  56 #ifdef MPI
  57       do iparm=1,nParmSet
  58         do ib=1,nT_h(iparm)
  59           do i=1,nR(ib,iparm)
  60             snk_p(i,ib,iparm)=0
  61           enddo
  62         enddo
  63       enddo
  64       do i=indstart(me1),indend(me1)
  65 #else
  66       do iparm=1,nParmSet
  67         do ib=1,nT_h(iparm)
  68           do i=1,nR(ib,iparm)
  69             snk(i,ib,iparm)=0
  70           enddo
  71         enddo
  72       enddo
  73       do i=1,ntot
  74 #endif
  75         read(ientout,rec=i,err=101)
  76      &    ((csingle(l,k),l=1,3),k=1,nres),
  77      &    ((csingle(l,k+nres),l=1,3),k=nnt,nct),
  78      &    nss,(ihpb(k),jhpb(k),k=1,nss),
  79      &    eini,efree,rmsdev,(q(j,iii+1),j=1,nQ),iR,ib,ipar
  80          if (indpdb.gt.0) then
  81            do k=1,nres
  82              do l=1,3
  83                c(l,k)=csingle(l,k)
  84              enddo
  85            enddo
  86            do k=nnt,nct
  87              do l=1,3
  88                c(l,k+nres)=csingle(l,k+nres)
  89              enddo
  90            enddo
  91            anatemp= 1.0d0/(beta_h(ib,ipar)*1.987D-3)
  92            q(nQ+1,iii+1)=rmsnat(iii+1)
  93          endif
  94          q(nQ+2,iii+1)=gyrate(iii+1)
  95 c        fT=T0*beta_h(ib,ipar)*1.987D-3
  96 c        ft=2.0d0/(1.0d0+1.0d0/(T0*beta_h(ib,ipar)*1.987D-3))
  97         if (rescale_mode.eq.1) then
  98           quot=1.0d0/(T0*beta_h(ib,ipar)*1.987D-3)
  99 #if defined(FUNCTH)
 100           tt = 1.0d0/(beta_h(ib,ipar)*1.987D-3)
 101           ft(6)=(320.0+80.0*dtanh((tt-320.0)/80.0))/320.0
 102 #elif defined(FUNCT)
 103           ft(6)=quot
 104 #else
 105           ft(6)=1.0d0
 106 #endif
 107           quotl=1.0d0
 108           kfacl=1.0d0
 109           do l=1,5
 110             quotl=quotl*quot
 111             kfacl=kfacl*kfac
 112             fT(l)=kfacl/(kfacl-1.0d0+quotl)
 113           enddo
 114         else if (rescale_mode.eq.2) then
 115           quot=1.0d0/(T0*beta_h(ib,ipar)*1.987D-3)
 116 #if defined(FUNCTH)
 117           tt = 1.0d0/(beta_h(ib,ipar)*1.987D-3)
 118           ft(6)=(320.0+80.0*dtanh((tt-320.0)/80.0))/320.0
 119 #elif defined(FUNCT)
 120           ft(6)=quot
 121 #else
 122           ft(6)=1.0d0
 123 #endif
 124           quotl=1.0d0
 125           do l=1,5
 126             quotl=quotl*quot
 127             fT(l)=1.12692801104297249644d0/
 128      &         dlog(dexp(quotl)+dexp(-quotl))
 129           enddo
 130         else if (rescale_mode.eq.0) then
 131           do l=1,5
 132             fT(l)=1.0d0
 133           enddo
 134         else
 135           write (iout,*) "Error in ECECALC: wrong RESCALE_MODE",
 136      &     rescale_mode
 137           call flush(iout)
 138           return1
 139         endif
 140
 141 c        write (iout,*) "T",1.0d0/(beta_h(ib,ipar)*1.987D-3)," T0",T0,
 142 c     &   " kfac",kfac,"quot",quot," fT",fT
 143         do j=1,2*nres
 144           do k=1,3
 145             c(k,j)=csingle(k,j)
 146           enddo
 147         enddo
 148         call int_from_cart1(.false.)
 149         ii=ii+1
 150         do iparm=1,nparmset
 151
 152         call restore_parm(iparm)
 153 #ifdef DEBUG
 154             write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
 155      &      wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
 156      &      wtor_d,wsccor,wbond
 157 #endif
 158         call etotal(energia(0),fT)
 159 #ifdef DEBUG
 160         write (iout,*) "Conformation",i
 161           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 162           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 163         call enerprint(energia(0),fT)
 164         write (iout,'(2i5,21f8.2)') i,iparm,(energia(k),k=1,21)
 165         write (iout,*) "ftors",ftors
 166         call briefout(i,energia(0))
 167         temp=1.0d0/(beta_h(ib,ipar)*1.987D-3)
 168         write (iout,*) "temp", temp
 169         call pdbout(i,temp,energia(0),energia(0),0.0d0,0.0d0)
 170 #endif
 171         if (energia(0).ge.1.0d20) then
 172           write (iout,*) "NaNs detected in some of the energy",
 173      &     " components for conformation",ii+1
 174           write (iout,*) "The Cartesian geometry is:"
 175           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 176           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 177           write (iout,*) "The internal geometry is:"
 178 c          call intout
 179 c        call pdbout(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
 180           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 181           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 182           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 183           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 184           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 185           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 186           write (iout,*) "The components of the energy are:"
 187           call enerprint(energia(0),fT)
 188           write (iout,*)
 189      &      "This conformation WILL NOT be added to the database."
 190           call flush(iout)
 191           goto 121
 192         else
 193 #ifdef DEBUG
 194           if (ipar.eq.iparm) write (iout,*) i,iparm,
 195      &      1.0d0/(beta_h(ib,ipar)*1.987D-3),eini,energia(0)
 196 #endif
 197           if (ipar.eq.iparm .and. einicheck.gt.0 .and.
 198      &      dabs(eini-energia(0)).gt.tole) then
 199             if (errmsg_count.le.maxerrmsg_count) then
 200               write (iout,'(2a,2e15.5,a,2i8,a,f8.1)')
 201      &         "Warning: energy differs remarkably from ",
 202      &         " the value read in: ",energia(0),eini," point",
 203      &         iii+1,indstart(me1)+iii," T",
 204      &         1.0d0/(1.987D-3*beta_h(ib,ipar))
 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       character*64 nazwa
 351       character*80 bxname,cxname
 352       character*64 bprotfile_temp
 353       character*3 liczba,licz
 354       character*2 licz2
 355       integer i,itj,ii,iii,j,k,l
 356       integer ixdrf,iret
 357       integer iscor,islice
 358       double precision rmsdev,efree,eini
 359       real*4 csingle(3,maxres2)
 360       double precision energ
 361       integer ilen,iroof
 362       external ilen,iroof
 363       integer ir,ib,iparm
 364       write (licz2,'(bz,i2.2)') islice
 365       call opentmp(islice,ientout,bprotfile_temp)
 366       write (iout,*) "bprotfile_temp ",bprotfile_temp
 367       call flush(iout)
 368       if (.not.bxfile .and. .not. cxfile .and. indpdb.eq.0
 369      &   .and. ensembles.eq.0) then
 370         close(ientout,status="delete")
 371         return
 372       endif
 373 #ifdef MPI
 374       write (liczba,'(bz,i3.3)') me
 375       if (bxfile .or. cxfile .or. ensembles.gt.0) then
 376         if (.not.separate_parset) then
 377           bxname = prefix(:ilen(prefix))//liczba//".bx"
 378         else
 379           write (licz,'(bz,i3.3)') myparm
 380           bxname = prefix(:ilen(prefix))//liczba//"_par"//licz//".bx"
 381         endif
 382         open (ientin,file=bxname,status="unknown",
 383      &    form="unformatted",access="direct",recl=lenrec1)
 384       endif
 385 #else
 386       if (bxfile .or. cxfile .or. ensembles.gt.0) then
 387         if (nslice.eq.1) then
 388           bxname = prefix(:ilen(prefix))//".bx"
 389         else
 390           bxname = prefix(:ilen(prefix))//
 391      &           "_slice_"//licz2//".bx"
 392         endif
 393         open (ientin,file=bxname,status="unknown",
 394      &    form="unformatted",access="direct",recl=lenrec1)
 395         write (iout,*) "Calculating energies; writing geometry",
 396      & " and energy components to ",bxname(:ilen(bxname))
 397       endif
 398 #if (defined(AIX) && !defined(JUBL))
 399         call xdrfopen_(ixdrf,cxname, "w", iret)
 400 #else
 401         call xdrfopen(ixdrf,cxname, "w", iret)
 402 #endif
 403         if (iret.eq.0) then
 404           write (iout,*) "Error opening cxfile ",cxname(:ilen(cxname))
 405           cxfile=.fale.
 406         endif
 407       endif
 408 #endif
 409       if (indpdb.gt.0) then
 410         if (nslice.eq.1) then
 411 #ifdef MPI
 412          if (.not.separate_parset) then
 413            statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))
 414      &       //liczba//'.stat'
 415          else
 416            write (licz,'(bz,i3.3)') myparm
 417            statname=prefix(:ilen(prefix))//'_par'//licz//'_'//
 418      &      pot(:ilen(pot))//liczba//'.stat'
 419          endif
 420
 421 #else
 422           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))//'.stat'
 423 #endif
 424         else
 425 #ifdef MPI
 426          if (.not.separate_parset) then
 427           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))//
 428      &      "_slice_"//licz2//liczba//'.stat'
 429          else
 430           write (licz,'(bz,i3.3)') myparm
 431           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))//
 432      &      '_par'//licz//"_slice_"//licz2//liczba//'.stat'
 433          endif
 434 #else
 435           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))
 436      &      //"_slice_"//licz2//'.stat'
 437 #endif
 438         endif
 439         open(istat,file=statname,status="unknown")
 440       endif
 441
 442 #ifdef MPI
 443       do i=1,scount(me)
 444 #else
 445       do i=1,ntot(islice)
 446 #endif
 447         read(ientout,rec=i,err=101)
 448      &    ((csingle(l,k),l=1,3),k=1,nres),
 449      &    ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 450      &    nss,(ihpb(k),jhpb(k),k=1,nss),
 451      &    eini,efree,rmsdev,(q(k,i),k=1,nQ),iR,ib,iparm
 452 c        write (iout,*) iR,ib,iparm,eini,efree
 453         do j=1,2*nres
 454           do k=1,3
 455             c(k,j)=csingle(k,j)
 456           enddo
 457         enddo
 458         call int_from_cart1(.false.)
 459         iscore=0
 460 c        write (iout,*) "Calling conf_compar",i
 461 c        call flush(iout)
 462         if (indpdb.gt.0) then
 463           call conf_compar(i,.false.,.true.)
 464         endif
 465 c        write (iout,*) "Exit conf_compar",i
 466 c        call flush(iout)
 467         if (bxfile .or.cxfile .or. ensembles.gt.0) write (ientin,rec=i)
 468      &    ((csingle(l,k),l=1,3),k=1,nres),
 469      &    ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 470      &    nss,(ihpb(k),jhpb(k),k=1,nss),
 471 c     &    potE(i,iparm),-entfac(i),rms_nat,iscore
 472      &    potE(i,nparmset),-entfac(i),rms_nat,iscore
 473 c        write (iout,'(2i5,3e15.5)') i,me,potE(i,iparm),-entfac(i)
 474 #ifndef MPI
 475         if (cxfile) call cxwrite(ixdrf,csingle,potE(i,nparmset),
 476      &    -entfac(i),rms_nat,iscore)
 477 #endif
 478       enddo
 479       close(ientout,status="delete")
 480       close(istat)
 481       if (bxfile .or. cxfile .or. ensembles.gt.0) close(ientin)
 482 #ifdef MPI
 483       call MPI_Barrier(WHAM_COMM,IERROR)
 484       if (me.ne.Master .or. .not.bxfile .and. .not. cxfile
 485      &   .and. ensembles.eq.0) return
 486       write (iout,*)
 487       if (bxfile .or. ensembles.gt.0) then
 488         if (nslice.eq.1) then
 489           if (.not.separate_parset) then
 490             bxname = prefix(:ilen(prefix))//".bx"
 491           else
 492             write (licz,'(bz,i3.3)') myparm
 493             bxname = prefix(:ilen(prefix))//"_par"//licz//".bx"
 494           endif
 495         else
 496           if (.not.separate_parset) then
 497             bxname = prefix(:ilen(prefix))//"_slice_"//licz2//".bx"
 498           else
 499             write (licz,'(bz,i3.3)') myparm
 500             bxname = prefix(:ilen(prefix))//"par_"//licz//
 501      &        "_slice_"//licz2//".bx"
 502           endif
 503         endif
 504         open (ientout,file=bxname,status="unknown",
 505      &      form="unformatted",access="direct",recl=lenrec1)
 506         write (iout,*) "Master is creating binary database ",
 507      &   bxname(:ilen(bxname))
 508       endif
 509       if (cxfile) then
 510         if (nslice.eq.1) then
 511           if (.not.separate_parset) then
 512             cxname = prefix(:ilen(prefix))//".cx"
 513           else
 514             cxname = prefix(:ilen(prefix))//"_par"//licz//".cx"
 515           endif
 516         else
 517           if (.not.separate_parset) then
 518             cxname = prefix(:ilen(prefix))//
 519      &             "_slice_"//licz2//".cx"
 520           else
 521             cxname = prefix(:ilen(prefix))//"_par"//licz//
 522      &             "_slice_"//licz2//".cx"
 523           endif
 524         endif
 525 #if (defined(AIX) && !defined(JUBL))
 526         call xdrfopen_(ixdrf,cxname, "w", iret)
 527 #else
 528         call xdrfopen(ixdrf,cxname, "w", iret)
 529 #endif
 530         if (iret.eq.0) then
 531           write (iout,*) "Error opening cxfile ",cxname(:ilen(cxname))
 532           cxfile=.false.
 533         endif
 534       endif
 535       do j=0,nprocs-1
 536         write (liczba,'(bz,i3.3)') j
 537         if (separate_parset) then
 538           write (licz,'(bz,i3.3)') myparm
 539           bxname = prefix(:ilen(prefix))//liczba//"_par"//licz//".bx"
 540         else
 541           bxname = prefix(:ilen(prefix))//liczba//".bx"
 542         endif
 543         open (ientin,file=bxname,status="unknown",
 544      &    form="unformatted",access="direct",recl=lenrec1)
 545         write (iout,*) "Master is reading conformations from ",
 546      &   bxname(:ilen(bxname))
 547         iii = 0
 548 c        write (iout,*) "j",j," indstart",indstart(j)," indend",indend(j)
 549 c        call flush(iout)
 550         do i=indstart(j),indend(j)
 551           iii = iii+1
 552           read(ientin,rec=iii,err=101)
 553      &      ((csingle(l,k),l=1,3),k=1,nres),
 554      &      ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 555      &      nss,(ihpb(k),jhpb(k),k=1,nss),
 556      &      eini,efree,rmsdev,iscor
 557           if (bxfile .or. ensembles.gt.0) then
 558             write (ientout,rec=i)
 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           endif
 564           if(cxfile)call cxwrite(ixdrf,csingle,eini,efree,rmsdev,iscor)
 565 #ifdef DEBUG
 566           do k=1,2*nres
 567             do l=1,3
 568               c(l,k)=csingle(l,k)
 569             enddo
 570           enddo
 571           call int_from_cart1(.false.)
 572           write (iout,'(2i5,3e15.5)') i,iii,eini,efree
 573           write (iout,*) "The Cartesian geometry is:"
 574           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 575           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 576           write (iout,*) "The internal geometry is:"
 577           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 578           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 579           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 580           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 581           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 582           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 583           write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
 584           write (iout,'(f10.5,i5)') rmsdev,iscor
 585 #endif
 586         enddo ! i
 587         write (iout,*) iii," conformations (from",indstart(j)," to",
 588      &   indend(j),") read from ",
 589      &   bxname(:ilen(bxname))
 590         close (ientin,status="delete")
 591       enddo ! j
 592       if (bxfile .or. cxfile .or. ensembles.gt.0) close (ientout)
 593 #if (defined(AIX) && !defined(JUBL))
 594       if (cxfile) call xdrfclose_(ixdrf,cxname,iret)
 595 #else
 596       if (cxfile) call xdrfclose(ixdrf,cxname,iret)
 597 #endif
 598 #endif
 599       return
 600   101 write (iout,*) "Error in scratchfile."
 601       call flush(iout)
 602       return1
 603       end
 604 c-------------------------------------------------------------------------------
 605       subroutine cxwrite(ixdrf,csingle,eini,efree,rmsdev,iscor)
 606       implicit none
 607       include "DIMENSIONS"
 608       include "DIMENSIONS.ZSCOPT"
 609       include "DIMENSIONS.FREE"
 610       include "DIMENSIONS.COMPAR"
 611 #ifdef MPI
 612       include "mpif.h"
 613       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 614       include "COMMON.MPI"
 615 #endif
 616       include "COMMON.CONTROL"
 617       include "COMMON.CHAIN"
 618       include "COMMON.IOUNITS"
 619       include "COMMON.PROTFILES"
 620       include "COMMON.NAMES"
 621       include "COMMON.VAR"
 622       include "COMMON.SBRIDGE"
 623       include "COMMON.GEO"
 624       include "COMMON.FFIELD"
 625       include "COMMON.ENEPS"
 626       include "COMMON.LOCAL"
 627       include "COMMON.WEIGHTS"
 628       include "COMMON.INTERACT"
 629       include "COMMON.FREE"
 630       include "COMMON.ENERGIES"
 631       include "COMMON.COMPAR"
 632       include "COMMON.PROT"
 633       integer i,j,itmp,iscor,iret,ixdrf
 634       double precision rmsdev,efree,eini
 635       real*4 csingle(3,maxres2),xoord(3,maxres2+2)
 636       real*4 prec
 637
 638 c      write (iout,*) "cxwrite"
 639 c      call flush(iout)
 640       prec=10000.0
 641       do i=1,nres
 642        do j=1,3
 643         xoord(j,i)=csingle(j,i)
 644        enddo
 645       enddo
 646       do i=nnt,nct
 647        do j=1,3
 648         xoord(j,nres+i-nnt+1)=csingle(j,i+nres)
 649        enddo
 650       enddo
 651
 652       itmp=nres+nct-nnt+1
 653
 654 c      write (iout,*) "itmp",itmp
 655 c      call flush(iout)
 656 #if (defined(AIX) && !defined(JUBL))
 657       call xdrf3dfcoord_(ixdrf, xoord, itmp, prec, iret)
 658
 659 c      write (iout,*) "xdrf3dfcoord"
 660 c      call flush(iout)
 661       call xdrfint_(ixdrf, nss, iret)
 662       do j=1,nss
 663         call xdrfint_(ixdrf, ihpb(j), iret)
 664         call xdrfint_(ixdrf, jhpb(j), iret)
 665       enddo
 666       call xdrffloat_(ixdrf,real(eini),iret)
 667       call xdrffloat_(ixdrf,real(efree),iret)
 668       call xdrffloat_(ixdrf,real(rmsdev),iret)
 669       call xdrfint_(ixdrf,iscor,iret)
 670 #else
 671       call xdrf3dfcoord(ixdrf, xoord, itmp, prec, iret)
 672
 673       call xdrfint(ixdrf, nss, iret)
 674       do j=1,nss
 675         call xdrfint(ixdrf, ihpb(j), iret)
 676         call xdrfint(ixdrf, jhpb(j), iret)
 677       enddo
 678       call xdrffloat(ixdrf,real(eini),iret)
 679       call xdrffloat(ixdrf,real(efree),iret)
 680       call xdrffloat(ixdrf,real(rmsdev),iret)
 681       call xdrfint(ixdrf,iscor,iret)
 682 #endif
 683
 684       return
 685       end
 686 c------------------------------------------------------------------------------
 687       logical function conf_check(ii,iprint)
 688       implicit none
 689       include "DIMENSIONS"
 690       include "DIMENSIONS.ZSCOPT"
 691       include "DIMENSIONS.FREE"
 692 #ifdef MPI
 693       include "mpif.h"
 694       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 695       include "COMMON.MPI"
 696 #endif
 697       include "COMMON.CHAIN"
 698       include "COMMON.IOUNITS"
 699       include "COMMON.PROTFILES"
 700       include "COMMON.NAMES"
 701       include "COMMON.VAR"
 702       include "COMMON.SBRIDGE"
 703       include "COMMON.GEO"
 704       include "COMMON.FFIELD"
 705       include "COMMON.ENEPS"
 706       include "COMMON.LOCAL"
 707       include "COMMON.WEIGHTS"
 708       include "COMMON.INTERACT"
 709       include "COMMON.FREE"
 710       include "COMMON.ENERGIES"
 711       include "COMMON.CONTROL"
 712       include "COMMON.TORCNSTR"
 713       integer j,k,l,ii,itj,iprint
 714       if (.not.check_conf) then
 715         conf_check=.true.
 716         return
 717       endif
 718       call int_from_cart1(.false.)
 719       do j=nnt+1,nct
 720         if (itype(j-1).ne.21 .and. itype(j).ne.21 .and.
 721      &    (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0)) then
 722           if (iprint.gt.0)
 723      &    write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),
 724      &      " for conformation",ii
 725           if (iprint.gt.1) then
 726             write (iout,*) "The Cartesian geometry is:"
 727             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 728             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 729             write (iout,*) "The internal geometry is:"
 730             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 731             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 732             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 733             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 734             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 735             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 736           endif
 737           if (iprint.gt.0) write (iout,*)
 738      &      "This conformation WILL NOT be added to the database."
 739           conf_check=.false.
 740           return
 741         endif
 742       enddo
 743       do j=nnt,nct
 744         itj=itype(j)
 745         if (itype(j).ne.10 .and.itype(j).ne.21 .and.
 746      &     (vbld(nres+j)-dsc(itj)).gt.2.0d0) then
 747           if (iprint.gt.0)
 748      &    write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j),
 749      &     " for conformation",ii
 750           if (iprint.gt.1) then
 751             write (iout,*) "The Cartesian geometry is:"
 752             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 753             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 754             write (iout,*) "The internal geometry is:"
 755             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 756             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 757             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 758             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 759             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 760             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 761           endif
 762           if (iprint.gt.0) write (iout,*)
 763      &      "This conformation WILL NOT be added to the database."
 764           conf_check=.false.
 765           return
 766         endif
 767       enddo
 768       do j=3,nres
 769         if (theta(j).le.0.0d0) then
 770           if (iprint.gt.0)
 771      &    write (iout,*) "Zero theta angle(s) in 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         if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
 790       enddo
 791       conf_check=.true.
 792 c      write (iout,*) "conf_check passed",ii
 793       return
 794       end