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