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