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