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