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,scount_buff(0:99)
  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(1)",ftors(1)
 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_buff(me)=iii
 284       write (iout,*) "Me",me," scount_buff",scount_buff(me)
 285       call flush(iout)
 286 c  Master gathers updated numbers of conformations written by all procs.
 287 c      call MPI_AllGather(MPI_IN_PLACE,1,MPI_DATATYPE_NULL,scount(0),1,
 288 c     &  MPI_INTEGER, WHAM_COMM, IERROR)
 289       call MPI_AllGather( scount_buff(me), 1, MPI_INTEGER, scount(0), 1,
 290      &  MPI_INTEGER, WHAM_COMM, IERROR)
 291
 292       indstart(0)=1
 293       indend(0)=scount(0)
 294       do i=1, Nprocs-1
 295         indstart(i)=indend(i-1)+1
 296         indend(i)=indstart(i)+scount(i)-1
 297       enddo
 298       write (iout,*)
 299       write (iout,*) "Revised conformation counts"
 300       do i=0,nprocs1-1
 301         write (iout,'(a,i5,a,i7,a,i7,a,i7)')
 302      &    "Processor",i," indstart",indstart(i),
 303      &    " indend",indend(i)," count",scount(i)
 304       enddo
 305       call flush(iout)
 306       call MPI_AllReduce(snk_p(1,1,1),snk(1,1,1,islice),
 307      &  MaxR*MaxT_h*nParmSet,
 308      &  MPI_INTEGER,MPI_SUM,WHAM_COMM,IERROR)
 309 #endif
 310       stot(islice)=0
 311       do iparm=1,nParmSet
 312         do ib=1,nT_h(iparm)
 313           do i=1,nR(ib,iparm)
 314             stot(islice)=stot(islice)+snk(i,ib,iparm,islice)
 315           enddo
 316         enddo
 317       enddo
 318       write (iout,*) "Revised SNK"
 319       do iparm=1,nParmSet
 320         do ib=1,nT_h(iparm)
 321           write (iout,'("Param",i3," Temp",f6.1,3x,32i8)')
 322      &     iparm,1.0d0/(1.987D-3*beta_h(ib,iparm)),
 323      &     (snk(i,ib,iparm,islice),i=1,nR(ib,iparm))
 324           write (iout,*) "snk_p",(snk_p(i,ib,iparm),i=1,nR(ib,iparm))
 325         enddo
 326       enddo
 327       write (iout,'("Total",i10)') stot(islice)
 328       call flush(iout)
 329       return
 330   101 write (iout,*) "Error in scratchfile."
 331       call flush(iout)
 332       return1
 333       end
 334 c------------------------------------------------------------------------------
 335       subroutine write_dbase(islice,*)
 336       implicit none
 337       include "DIMENSIONS"
 338       include "DIMENSIONS.ZSCOPT"
 339       include "DIMENSIONS.FREE"
 340       include "DIMENSIONS.COMPAR"
 341 #ifdef MPI
 342       include "mpif.h"
 343       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 344       include "COMMON.MPI"
 345 #endif
 346       include "COMMON.CONTROL"
 347       include "COMMON.CHAIN"
 348       include "COMMON.IOUNITS"
 349       include "COMMON.PROTFILES"
 350       include "COMMON.NAMES"
 351       include "COMMON.VAR"
 352       include "COMMON.SBRIDGE"
 353       include "COMMON.GEO"
 354       include "COMMON.FFIELD"
 355       include "COMMON.ENEPS"
 356       include "COMMON.LOCAL"
 357       include "COMMON.WEIGHTS"
 358       include "COMMON.INTERACT"
 359       include "COMMON.FREE"
 360       include "COMMON.ENERGIES"
 361       include "COMMON.COMPAR"
 362       include "COMMON.PROT"
 363       include "COMMON.CONTACTS1"
 364       character*64 nazwa
 365       character*80 bxname,cxname
 366       character*64 bprotfile_temp
 367       character*3 liczba,licz
 368       character*2 licz2
 369       integer i,itj,ii,iii,j,k,l
 370       integer ixdrf,iret
 371       integer iscor,islice
 372       double precision rmsdev,efree,eini
 373       real*4 csingle(3,maxres2)
 374       double precision energ
 375       integer ilen,iroof
 376       external ilen,iroof
 377       integer ir,ib,iparm, scount_buff(0:99)
 378       integer isecstr(maxres)
 379       write (licz2,'(bz,i2.2)') islice
 380       call opentmp(islice,ientout,bprotfile_temp)
 381       write (iout,*) "bprotfile_temp ",bprotfile_temp
 382       call flush(iout)
 383       if (.not.bxfile .and. .not. cxfile .and. indpdb.eq.0
 384      &   .and. ensembles.eq.0) then
 385         close(ientout,status="delete")
 386         return
 387       endif
 388 #ifdef MPI
 389       write (liczba,'(bz,i3.3)') me
 390       if (bxfile .or. cxfile .or. ensembles.gt.0) then
 391         if (.not.separate_parset) then
 392           bxname = prefix(:ilen(prefix))//liczba//".bx"
 393         else
 394           write (licz,'(bz,i3.3)') myparm
 395           bxname = prefix(:ilen(prefix))//liczba//"_par"//licz//".bx"
 396         endif
 397         open (ientin,file=bxname,status="unknown",
 398      &    form="unformatted",access="direct",recl=lenrec1)
 399       endif
 400 #else
 401       if (bxfile .or. cxfile .or. ensembles.gt.0) then
 402         if (nslice.eq.1) then
 403           bxname = prefix(:ilen(prefix))//".bx"
 404         else
 405           bxname = prefix(:ilen(prefix))//
 406      &           "_slice_"//licz2//".bx"
 407         endif
 408         open (ientin,file=bxname,status="unknown",
 409      &    form="unformatted",access="direct",recl=lenrec1)
 410         write (iout,*) "Calculating energies; writing geometry",
 411      & " and energy components to ",bxname(:ilen(bxname))
 412       endif
 413 #if (defined(AIX) && !defined(JUBL))
 414         call xdrfopen_(ixdrf,cxname, "w", iret)
 415 #else
 416         call xdrfopen(ixdrf,cxname, "w", iret)
 417 #endif
 418         if (iret.eq.0) then
 419           write (iout,*) "Error opening cxfile ",cxname(:ilen(cxname))
 420           cxfile=.fale.
 421         endif
 422       endif
 423 #endif
 424       if (indpdb.gt.0) then
 425         if (nslice.eq.1) then
 426 #ifdef MPI
 427          if (.not.separate_parset) then
 428            statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))
 429      &       //liczba//'.stat'
 430          else
 431            write (licz,'(bz,i3.3)') myparm
 432            statname=prefix(:ilen(prefix))//'_par'//licz//'_'//
 433      &      pot(:ilen(pot))//liczba//'.stat'
 434          endif
 435
 436 #else
 437           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))//'.stat'
 438 #endif
 439         else
 440 #ifdef MPI
 441          if (.not.separate_parset) then
 442           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))//
 443      &      "_slice_"//licz2//liczba//'.stat'
 444          else
 445           write (licz,'(bz,i3.3)') myparm
 446           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))//
 447      &      '_par'//licz//"_slice_"//licz2//liczba//'.stat'
 448          endif
 449 #else
 450           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))
 451      &      //"_slice_"//licz2//'.stat'
 452 #endif
 453         endif
 454         open(istat,file=statname,status="unknown")
 455       endif
 456
 457 #ifdef MPI
 458       do i=1,scount(me)
 459 #else
 460       do i=1,ntot(islice)
 461 #endif
 462         read(ientout,rec=i,err=101)
 463      &    ((csingle(l,k),l=1,3),k=1,nres),
 464      &    ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 465      &    nss,(ihpb(k),jhpb(k),k=1,nss),
 466      &    eini,efree,rmsdev,(q(k,i),k=1,nQ),iR,ib,iparm
 467 c        write (iout,*) iR,ib,iparm,eini,efree
 468         do j=1,2*nres
 469           do k=1,3
 470             c(k,j)=csingle(k,j)
 471           enddo
 472         enddo
 473         call int_from_cart1(.false.)
 474         iscore=0
 475 c        write (iout,*) "Calling conf_compar",i
 476 c        call flush(iout)
 477          anatemp= 1.0d0/(beta_h(ib,iparm)*1.987D-3)
 478         if (indpdb.gt.0) then
 479           call conf_compar(i,.false.,.true.)
 480 c        else
 481 c            call elecont(.false.,ncont,icont,nnt,nct)
 482 c            call secondary2(.false.,.false.,ncont,icont,isecstr)
 483         endif
 484 c        write (iout,*) "Exit conf_compar",i
 485 c        call flush(iout)
 486         if (bxfile .or.cxfile .or. ensembles.gt.0) write (ientin,rec=i)
 487      &    ((csingle(l,k),l=1,3),k=1,nres),
 488      &    ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 489      &    nss,(ihpb(k),jhpb(k),k=1,nss),
 490 c     &    potE(i,iparm),-entfac(i),rms_nat,iscore
 491      &    potE(i,nparmset),-entfac(i),rms_nat,iscore
 492 c        write (iout,'(2i5,3e15.5)') i,me,potE(i,iparm),-entfac(i)
 493 #ifndef MPI
 494         if (cxfile) call cxwrite(ixdrf,csingle,potE(i,nparmset),
 495      &    -entfac(i),rms_nat,iscore)
 496 #endif
 497       enddo
 498       close(ientout,status="delete")
 499       close(istat)
 500       if (bxfile .or. cxfile .or. ensembles.gt.0) close(ientin)
 501 #ifdef MPI
 502       call MPI_Barrier(WHAM_COMM,IERROR)
 503       if (me.ne.Master .or. .not.bxfile .and. .not. cxfile
 504      &   .and. ensembles.eq.0) return
 505       write (iout,*)
 506       if (bxfile .or. ensembles.gt.0) then
 507         if (nslice.eq.1) then
 508           if (.not.separate_parset) then
 509             bxname = prefix(:ilen(prefix))//".bx"
 510           else
 511             write (licz,'(bz,i3.3)') myparm
 512             bxname = prefix(:ilen(prefix))//"_par"//licz//".bx"
 513           endif
 514         else
 515           if (.not.separate_parset) then
 516             bxname = prefix(:ilen(prefix))//"_slice_"//licz2//".bx"
 517           else
 518             write (licz,'(bz,i3.3)') myparm
 519             bxname = prefix(:ilen(prefix))//"par_"//licz//
 520      &        "_slice_"//licz2//".bx"
 521           endif
 522         endif
 523         open (ientout,file=bxname,status="unknown",
 524      &      form="unformatted",access="direct",recl=lenrec1)
 525         write (iout,*) "Master is creating binary database ",
 526      &   bxname(:ilen(bxname))
 527       endif
 528       if (cxfile) then
 529         if (nslice.eq.1) then
 530           if (.not.separate_parset) then
 531             cxname = prefix(:ilen(prefix))//".cx"
 532           else
 533             cxname = prefix(:ilen(prefix))//"_par"//licz//".cx"
 534           endif
 535         else
 536           if (.not.separate_parset) then
 537             cxname = prefix(:ilen(prefix))//
 538      &             "_slice_"//licz2//".cx"
 539           else
 540             cxname = prefix(:ilen(prefix))//"_par"//licz//
 541      &             "_slice_"//licz2//".cx"
 542           endif
 543         endif
 544 #if (defined(AIX) && !defined(JUBL))
 545         call xdrfopen_(ixdrf,cxname, "w", iret)
 546 #else
 547         call xdrfopen(ixdrf,cxname, "w", iret)
 548 #endif
 549         if (iret.eq.0) then
 550           write (iout,*) "Error opening cxfile ",cxname(:ilen(cxname))
 551           cxfile=.false.
 552         endif
 553       endif
 554       do j=0,nprocs-1
 555         write (liczba,'(bz,i3.3)') j
 556         if (separate_parset) then
 557           write (licz,'(bz,i3.3)') myparm
 558           bxname = prefix(:ilen(prefix))//liczba//"_par"//licz//".bx"
 559         else
 560           bxname = prefix(:ilen(prefix))//liczba//".bx"
 561         endif
 562         open (ientin,file=bxname,status="unknown",
 563      &    form="unformatted",access="direct",recl=lenrec1)
 564         write (iout,*) "Master is reading conformations from ",
 565      &   bxname(:ilen(bxname))
 566         iii = 0
 567 c        write (iout,*) "j",j," indstart",indstart(j)," indend",indend(j)
 568 c        call flush(iout)
 569         do i=indstart(j),indend(j)
 570           iii = iii+1
 571           read(ientin,rec=iii,err=101)
 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           if (bxfile .or. ensembles.gt.0) then
 577             write (ientout,rec=i)
 578      &        ((csingle(l,k),l=1,3),k=1,nres),
 579      &        ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 580      &        nss,(ihpb(k),jhpb(k),k=1,nss),
 581      &        eini,efree,rmsdev,iscor
 582           endif
 583           if(cxfile)call cxwrite(ixdrf,csingle,eini,efree,rmsdev,iscor)
 584 #ifdef DEBUG
 585           do k=1,2*nres
 586             do l=1,3
 587               c(l,k)=csingle(l,k)
 588             enddo
 589           enddo
 590           call int_from_cart1(.false.)
 591           write (iout,'(2i5,3e15.5)') i,iii,eini,efree
 592           write (iout,*) "The Cartesian geometry is:"
 593           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 594           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 595           write (iout,*) "The internal geometry is:"
 596           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 597           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 598           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 599           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 600           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 601           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 602           write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
 603           write (iout,'(f10.5,i5)') rmsdev,iscor
 604 #endif
 605         enddo ! i
 606         write (iout,*) iii," conformations (from",indstart(j)," to",
 607      &   indend(j),") read from ",
 608      &   bxname(:ilen(bxname))
 609         close (ientin,status="delete")
 610       enddo ! j
 611       if (bxfile .or. cxfile .or. ensembles.gt.0) close (ientout)
 612 #if (defined(AIX) && !defined(JUBL))
 613       if (cxfile) call xdrfclose_(ixdrf,cxname,iret)
 614 #else
 615       if (cxfile) call xdrfclose(ixdrf,cxname,iret)
 616 #endif
 617 #endif
 618       return
 619   101 write (iout,*) "Error in scratchfile."
 620       call flush(iout)
 621       return1
 622       end
 623 c-------------------------------------------------------------------------------
 624       subroutine cxwrite(ixdrf,csingle,eini,efree,rmsdev,iscor)
 625       implicit none
 626       include "DIMENSIONS"
 627       include "DIMENSIONS.ZSCOPT"
 628       include "DIMENSIONS.FREE"
 629       include "DIMENSIONS.COMPAR"
 630 #ifdef MPI
 631       include "mpif.h"
 632       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 633       include "COMMON.MPI"
 634 #endif
 635       include "COMMON.CONTROL"
 636       include "COMMON.CHAIN"
 637       include "COMMON.IOUNITS"
 638       include "COMMON.PROTFILES"
 639       include "COMMON.NAMES"
 640       include "COMMON.VAR"
 641       include "COMMON.SBRIDGE"
 642       include "COMMON.GEO"
 643       include "COMMON.FFIELD"
 644       include "COMMON.ENEPS"
 645       include "COMMON.LOCAL"
 646       include "COMMON.WEIGHTS"
 647       include "COMMON.INTERACT"
 648       include "COMMON.FREE"
 649       include "COMMON.ENERGIES"
 650       include "COMMON.COMPAR"
 651       include "COMMON.PROT"
 652       integer i,j,itmp,iscor,iret,ixdrf
 653       double precision rmsdev,efree,eini
 654       real*4 csingle(3,maxres2),xoord(3,maxres2+2)
 655       real*4 prec
 656
 657 c      write (iout,*) "cxwrite"
 658 c      call flush(iout)
 659       prec=10000.0
 660       do i=1,nres
 661        do j=1,3
 662         xoord(j,i)=csingle(j,i)
 663        enddo
 664       enddo
 665       do i=nnt,nct
 666        do j=1,3
 667         xoord(j,nres+i-nnt+1)=csingle(j,i+nres)
 668        enddo
 669       enddo
 670
 671       itmp=nres+nct-nnt+1
 672
 673 c      write (iout,*) "itmp",itmp
 674 c      call flush(iout)
 675 #if (defined(AIX) && !defined(JUBL))
 676       call xdrf3dfcoord_(ixdrf, xoord, itmp, prec, iret)
 677
 678 c      write (iout,*) "xdrf3dfcoord"
 679 c      call flush(iout)
 680       call xdrfint_(ixdrf, nss, iret)
 681       do j=1,nss
 682            if (dyn_ss) then
 683             call xdrfint(ixdrf, idssb(j)+nres, iret)
 684             call xdrfint(ixdrf, jdssb(j)+nres, iret)
 685            else
 686             call xdrfint_(ixdrf, ihpb(j), iret)
 687             call xdrfint_(ixdrf, jhpb(j), iret)
 688            endif
 689       enddo
 690       call xdrffloat_(ixdrf,real(eini),iret)
 691       call xdrffloat_(ixdrf,real(efree),iret)
 692       call xdrffloat_(ixdrf,real(rmsdev),iret)
 693       call xdrfint_(ixdrf,iscor,iret)
 694 #else
 695       call xdrf3dfcoord(ixdrf, xoord, itmp, prec, iret)
 696
 697       call xdrfint(ixdrf, nss, iret)
 698       do j=1,nss
 699            if (dyn_ss) then
 700             call xdrfint(ixdrf, idssb(j)+nres, iret)
 701             call xdrfint(ixdrf, jdssb(j)+nres, iret)
 702            else
 703             call xdrfint(ixdrf, ihpb(j), iret)
 704             call xdrfint(ixdrf, jhpb(j), iret)
 705            endif
 706       enddo
 707       call xdrffloat(ixdrf,real(eini),iret)
 708       call xdrffloat(ixdrf,real(efree),iret)
 709       call xdrffloat(ixdrf,real(rmsdev),iret)
 710       call xdrfint(ixdrf,iscor,iret)
 711 #endif
 712
 713       return
 714       end
 715 c------------------------------------------------------------------------------
 716       logical function conf_check(ii,iprint)
 717       implicit none
 718       include "DIMENSIONS"
 719       include "DIMENSIONS.ZSCOPT"
 720       include "DIMENSIONS.FREE"
 721 #ifdef MPI
 722       include "mpif.h"
 723       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 724       include "COMMON.MPI"
 725 #endif
 726       include "COMMON.CHAIN"
 727       include "COMMON.IOUNITS"
 728       include "COMMON.PROTFILES"
 729       include "COMMON.NAMES"
 730       include "COMMON.VAR"
 731       include "COMMON.SBRIDGE"
 732       include "COMMON.GEO"
 733       include "COMMON.FFIELD"
 734       include "COMMON.ENEPS"
 735       include "COMMON.LOCAL"
 736       include "COMMON.WEIGHTS"
 737       include "COMMON.INTERACT"
 738       include "COMMON.FREE"
 739       include "COMMON.ENERGIES"
 740       include "COMMON.CONTROL"
 741       include "COMMON.TORCNSTR"
 742       integer j,k,l,ii,itj,iprint
 743       if (.not.check_conf) then
 744         conf_check=.true.
 745         return
 746       endif
 747       call int_from_cart1(.false.)
 748       do j=nnt+1,nct
 749         if (wliptran.gt.0d0) then
 750         if (itype(j-1).ne.ntyp1 .and. itype(j).ne.ntyp1 .and.
 751      &    (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.3d0)) then
 752           if (iprint.gt.0)
 753      &    write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),
 754      &      " for conformation",ii,wliptran
 755           if (iprint.gt.1) then
 756             write (iout,*) "The Cartesian geometry is:"
 757             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 758             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 759             write (iout,*) "The internal geometry is:"
 760             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 761             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 762             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 763             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 764             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 765             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 766           endif
 767           if (iprint.gt.0) write (iout,*)
 768      &      "This conformation WILL NOT be added to the database."
 769           conf_check=.false.
 770           return
 771         endif
 772
 773         else
 774         if (itype(j-1).ne.ntyp1 .and. itype(j).ne.ntyp1 .and.
 775      &    (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0)) then
 776           if (iprint.gt.0)
 777      &    write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),
 778      &      " for conformation",ii,wliptran
 779           if (iprint.gt.1) then
 780             write (iout,*) "The Cartesian geometry is:"
 781             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 782             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 783             write (iout,*) "The internal geometry is:"
 784             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 785             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 786             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 787             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 788             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 789             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 790           endif
 791           if (iprint.gt.0) write (iout,*)
 792      &      "This conformation WILL NOT be added to the database."
 793           conf_check=.false.
 794           return
 795         endif
 796         endif
 797       enddo
 798       do j=nnt,nct
 799         itj=itype(j)
 800         if (itype(j).ne.10 .and.itype(j).ne.ntyp1 .and.
 801      &     (vbld(nres+j)-dsc(iabs(itj))).gt.2.0d0) then
 802           if (iprint.gt.0)
 803      &    write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j),
 804      &     " for conformation",ii
 805           if (iprint.gt.1) then
 806             write (iout,*) "The Cartesian geometry is:"
 807             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 808             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 809             write (iout,*) "The internal geometry is:"
 810             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 811             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 812             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 813             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 814             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 815             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 816           endif
 817           if (iprint.gt.0) write (iout,*)
 818      &      "This conformation WILL NOT be added to the database."
 819           conf_check=.false.
 820           return
 821         endif
 822       enddo
 823       do j=3,nres
 824         if (theta(j).le.0.0d0) then
 825           if (iprint.gt.0)
 826      &    write (iout,*) "Zero theta angle(s) in conformation",ii
 827           if (iprint.gt.1) then
 828             write (iout,*) "The Cartesian geometry is:"
 829             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 830             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 831             write (iout,*) "The internal geometry is:"
 832             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 833             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 834             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 835             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 836             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 837             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 838           endif
 839           if (iprint.gt.0) write (iout,*)
 840      &      "This conformation WILL NOT be added to the database."
 841           conf_check=.false.
 842           return
 843         endif
 844         if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
 845       enddo
 846       conf_check=.true.
 847 c      write (iout,*) "conf_check passed",ii
 848       return
 849       end