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