source/unres/src_MD-M-SAXS/readrtns_CSA.F

   1       subroutine readrtns
   2       implicit real*8 (a-h,o-z)
   3       include 'DIMENSIONS'
   4 #ifdef MPI
   5       include 'mpif.h'
   6 #endif
   7       include 'COMMON.SETUP'
   8       include 'COMMON.CONTROL'
   9       include 'COMMON.SBRIDGE'
  10       include 'COMMON.IOUNITS'
  11       include 'COMMON.SPLITELE'
  12       logical file_exist
  13 C Read job setup parameters
  14       call read_control
  15 C Read force-field parameters except weights
  16       call parmread
  17 C Read control parameters for energy minimzation if required
  18       if (minim) call read_minim
  19 C Read MCM control parameters if required
  20       if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
  21 C Read MD control parameters if reqjuired
  22       if (modecalc.eq.12) call read_MDpar
  23 C Read MREMD control parameters if required
  24       if (modecalc.eq.14) then
  25          call read_MDpar
  26          call read_REMDpar
  27       endif
  28 C Read MUCA control parameters if required
  29       if (lmuca) call read_muca
  30 C Read CSA control parameters if required (from fort.40 if exists
  31 C otherwise from general input file)
  32       if (modecalc.eq.8) then
  33        inquire (file="fort.40",exist=file_exist)
  34        if (.not.file_exist) call csaread
  35       endif
  36 cfmc      if (modecalc.eq.10) call mcmfread
  37 c Read energy-term weights and disulfide parameters
  38       call weightread
  39 C Read molecule information, molecule geometry, energy-term weights, and
  40 C restraints if requested
  41       call molread
  42 C Print restraint information
  43 #ifdef MPI
  44       if (.not. out1file .or. me.eq.king) then
  45 #endif
  46       if (nhpb.gt.nss)  then
  47         write (iout,'(/a,i5,a/4a5,2a8,3a10,2a5)')
  48      &  "The following",nhpb-nss,
  49      &  " distance restraints have been imposed:",
  50      &  "   Nr"," res1"," res2"," beta","   d1","   d2","    k","    V",
  51      &  "  score"," type"
  52         do i=nss+1,nhpb
  53           write (iout,'(4i5,2f8.2,3f10.5,2i5)')i-nss,ihpb(i),jhpb(i),
  54      &  ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
  55      &  irestr_type(i)
  56         enddo
  57       endif
  58       if (npeak.gt.0) then
  59         write (iout,'(/a,i5,a/4a5,2a8,3a10,2a5)')
  60      &  "The following",npeak,
  61      &  " NMR peak restraints have been imposed:",
  62      &  "   Nr"," res1"," res2"," beta","   d1","   d2","    k","    V",
  63      &  "  score"," type"," ipeak"
  64         do i=1,npeak
  65           do j=ipeak(1,i),ipeak(2,i)
  66             write (iout,'(5i5,2f8.2,2f10.5,i5)')i,j,ihpb_peak(j),
  67      &       jhpb_peak(j),ibecarb_peak(j),dhpb_peak(j),dhpb1_peak(j),
  68      &       forcon_peak(j),fordepth_peak(i),irestr_type_peak(j)
  69           enddo
  70         enddo
  71         write (iout,*) "The ipeak array"
  72         do i=1,npeak
  73           write (iout,'(3i5)' ) i,ipeak(1,i),ipeak(2,i)
  74         enddo
  75       endif
  76 #ifdef MPI
  77       endif
  78 #endif
  79 c      print *,"Processor",myrank," leaves READRTNS"
  80       return
  81       end
  82 C-------------------------------------------------------------------------------
  83       subroutine read_control
  84 C
  85 C Read contorl data
  86 C
  87       implicit real*8 (a-h,o-z)
  88       include 'DIMENSIONS'
  89 #ifdef MP
  90       include 'mpif.h'
  91       logical OKRandom, prng_restart
  92       real*8  r1
  93 #endif
  94       include 'COMMON.IOUNITS'
  95       include 'COMMON.TIME1'
  96       include 'COMMON.THREAD'
  97       include 'COMMON.SBRIDGE'
  98       include 'COMMON.CONTROL'
  99       include 'COMMON.MCM'
 100       include 'COMMON.MAP'
 101       include 'COMMON.HEADER'
 102       include 'COMMON.CSA'
 103       include 'COMMON.CHAIN'
 104       include 'COMMON.MUCA'
 105       include 'COMMON.MD'
 106       include 'COMMON.FFIELD'
 107       include 'COMMON.INTERACT'
 108       include 'COMMON.SETUP'
 109       include 'COMMON.SPLITELE'
 110       include 'COMMON.SHIELD'
 111       include 'COMMON.GEO'
 112       COMMON /MACHSW/ KDIAG,ICORFL,IXDR
 113       character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
 114       character*80 ucase
 115       character*320 controlcard
 116
 117       nglob_csa=0
 118       eglob_csa=1d99
 119       nmin_csa=0
 120       read (INP,'(a)') titel
 121       call card_concat(controlcard)
 122 c      out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
 123 c      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
 124       call reada(controlcard,'SEED',seed,0.0D0)
 125       call random_init(seed)
 126 C Set up the time limit (caution! The time must be input in minutes!)
 127       read_cart=index(controlcard,'READ_CART').gt.0
 128       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
 129 C this variable with_theta_constr is the variable which allow to read and execute the
 130 C constrains on theta angles WITH_THETA_CONSTR is the keyword
 131       with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
 132       write (iout,*) "with_theta_constr ",with_theta_constr
 133       call readi(controlcard,'NSAXS',nsaxs,0)
 134       call readi(controlcard,'SAXS_MODE',saxs_mode,0)
 135       call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0)
 136       call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0)
 137       write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
 138      &   SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff
 139       call readi(controlcard,'SYM',symetr,1)
 140       call reada(controlcard,'TIMLIM',timlim,2800.0D0) ! default 16 hours
 141       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
 142       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
 143       call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
 144       call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
 145       call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
 146       call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
 147       call reada(controlcard,'DRMS',drms,0.1D0)
 148       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 149        write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
 150        write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
 151        write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
 152        write (iout,'(a,f10.1)')'DRMS    = ',drms
 153        write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
 154        write (iout,'(a,f10.1)') 'Time limit (min):',timlim
 155       endif
 156       call readi(controlcard,'NZ_START',nz_start,0)
 157       call readi(controlcard,'NZ_END',nz_end,0)
 158 c      call readi(controlcard,'IZ_SC',iz_sc,0)
 159       timlim=60.0D0*timlim
 160       safety = 60.0d0*safety
 161       timem=timlim
 162       modecalc=0
 163       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 164       minim=(index(controlcard,'MINIMIZE').gt.0)
 165       dccart=(index(controlcard,'CART').gt.0)
 166       overlapsc=(index(controlcard,'OVERLAP').gt.0)
 167       overlapsc=.not.overlapsc
 168       searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
 169       searchsc=.not.searchsc
 170       sideadd=(index(controlcard,'SIDEADD').gt.0)
 171       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
 172       mremd_dec=(index(controlcard,'MREMD_DEC').gt.0)
 173       outpdb=(index(controlcard,'PDBOUT').gt.0)
 174       outmol2=(index(controlcard,'MOL2OUT').gt.0)
 175       pdbref=(index(controlcard,'PDBREF').gt.0)
 176       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
 177       indpdb=index(controlcard,'PDBSTART')
 178       extconf=(index(controlcard,'EXTCONF').gt.0)
 179       AFMlog=(index(controlcard,'AFM'))
 180       selfguide=(index(controlcard,'SELFGUIDE'))
 181 c      print *,'AFMlog',AFMlog,selfguide,"KUPA"
 182       call readi(controlcard,'TUBEMOD',tubelog,0)
 183 c      write (iout,*) TUBElog,"TUBEMODE"
 184       call readi(controlcard,'IPRINT',iprint,0)
 185 C SHIELD keyword sets if the shielding effect of side-chains is used
 186 C 0 denots no shielding is used all peptide are equally despite the
 187 C solvent accesible area
 188 C 1 the newly introduced function
 189 C 2 reseved for further possible developement
 190       call readi(controlcard,'SHIELD',shield_mode,0)
 191 C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 192         write(iout,*) "shield_mode",shield_mode
 193 C      endif
 194       call readi(controlcard,'TORMODE',tor_mode,0)
 195 C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 196         write(iout,*) "torsional and valence angle mode",tor_mode
 197       call readi(controlcard,'MAXGEN',maxgen,10000)
 198       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
 199       call readi(controlcard,"KDIAG",kdiag,0)
 200       call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
 201       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
 202      & write (iout,*) "RESCALE_MODE",rescale_mode
 203       split_ene=index(controlcard,'SPLIT_ENE').gt.0
 204       if (index(controlcard,'REGULAR').gt.0.0D0) then
 205         call reada(controlcard,'WEIDIS',weidis,0.1D0)
 206         modecalc=1
 207         refstr=.true.
 208       endif
 209       if (refstr) then
 210       if (index(controlcard,"CASC").gt.0) then
 211         iz_sc=1
 212       else if (index(controlcard,"CAONLY").gt.0) then
 213         iz_sc=0
 214       else if (index(controlcard,"SCONLY").gt.0) then
 215         iz_sc=2
 216       else
 217         write (iout,*) "ERROR! Must specify CASC CAONLY or SCONLY when",
 218      &    " specifying REFSTR, PDBREF or REGULAR."
 219         stop
 220       endif
 221       endif
 222       if (index(controlcard,'CHECKGRAD').gt.0) then
 223         modecalc=5
 224         if (index(controlcard,' CART').gt.0) then
 225           icheckgrad=1
 226         elseif (index(controlcard,'CARINT').gt.0) then
 227           icheckgrad=2
 228         else
 229           icheckgrad=3
 230         endif
 231         call reada(controlcard,'DELTA',aincr,1.0d-7)
 232 c        write (iout,*) "icheckgrad",icheckgrad
 233       elseif (index(controlcard,'THREAD').gt.0) then
 234         modecalc=2
 235         call readi(controlcard,'THREAD',nthread,0)
 236         if (nthread.gt.0) then
 237           call reada(controlcard,'WEIDIS',weidis,0.1D0)
 238         else
 239           if (fg_rank.eq.0)
 240      &    write (iout,'(a)')'A number has to follow the THREAD keyword.'
 241           stop 'Error termination in Read_Control.'
 242         endif
 243       else if (index(controlcard,'MCMA').gt.0) then
 244         modecalc=3
 245       else if (index(controlcard,'MCEE').gt.0) then
 246         modecalc=6
 247       else if (index(controlcard,'MULTCONF').gt.0) then
 248         modecalc=4
 249       else if (index(controlcard,'MAP').gt.0) then
 250         modecalc=7
 251         call readi(controlcard,'MAP',nmap,0)
 252       else if (index(controlcard,'CSA').gt.0) then
 253         modecalc=8
 254 crc      else if (index(controlcard,'ZSCORE').gt.0) then
 255 crc
 256 crc  ZSCORE is rm from UNRES, modecalc=9 is available
 257 crc
 258 crc        modecalc=9
 259 cfcm      else if (index(controlcard,'MCMF').gt.0) then
 260 cfmc        modecalc=10
 261       else if (index(controlcard,'SOFTREG').gt.0) then
 262         modecalc=11
 263       else if (index(controlcard,'CHECK_BOND').gt.0) then
 264         modecalc=-1
 265       else if (index(controlcard,'TEST').gt.0) then
 266         modecalc=-2
 267       else if (index(controlcard,'MD').gt.0) then
 268         modecalc=12
 269       else if (index(controlcard,'RE ').gt.0) then
 270         modecalc=14
 271       endif
 272
 273       lmuca=index(controlcard,'MUCA').gt.0
 274       call readi(controlcard,'MUCADYN',mucadyn,0)
 275       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
 276       if (lmuca .and. (me.eq.king .or. .not.out1file ))
 277      & then
 278        write (iout,*) 'MUCADYN=',mucadyn
 279        write (iout,*) 'MUCASMOOTH=',muca_smooth
 280       endif
 281
 282       iscode=index(controlcard,'ONE_LETTER')
 283       indphi=index(controlcard,'PHI')
 284       indback=index(controlcard,'BACK')
 285       iranconf=index(controlcard,'RAND_CONF')
 286       i2ndstr=index(controlcard,'USE_SEC_PRED')
 287       gradout=index(controlcard,'GRADOUT').gt.0
 288       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
 289 C  DISTCHAINMAX become obsolete for periodic boundry condition
 290       call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
 291 C Reading the dimensions of box in x,y,z coordinates
 292       call reada(controlcard,'BOXX',boxxsize,100.0d0)
 293       call reada(controlcard,'BOXY',boxysize,100.0d0)
 294       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
 295 c Cutoff range for interactions
 296       call reada(controlcard,"R_CUT",r_cut,15.0d0)
 297       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
 298       call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
 299       call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
 300       if (lipthick.gt.0.0d0) then
 301        bordliptop=(boxzsize+lipthick)/2.0
 302        bordlipbot=bordliptop-lipthick
 303 C      endif
 304       if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
 305      & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
 306       buflipbot=bordlipbot+lipbufthick
 307       bufliptop=bordliptop-lipbufthick
 308       if ((lipbufthick*2.0d0).gt.lipthick)
 309      &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
 310       endif
 311       write(iout,*) "bordliptop=",bordliptop
 312       write(iout,*) "bordlipbot=",bordlipbot
 313       write(iout,*) "bufliptop=",bufliptop
 314       write(iout,*) "buflipbot=",buflipbot
 315       write (iout,*) "SHIELD MODE",shield_mode
 316       if (TUBElog.gt.0) then
 317        call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
 318        call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
 319        call reada(controlcard,"RTUBE",tubeR0,0.0d0)
 320        call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
 321        call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
 322        call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
 323        buftubebot=bordtubebot+tubebufthick
 324        buftubetop=bordtubetop-tubebufthick
 325       endif
 326 c      if (shield_mode.gt.0) then
 327 c      pi=3.141592d0
 328 C VSolvSphere the volume of solving sphere
 329 C      print *,pi,"pi"
 330 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
 331 C there will be no distinction between proline peptide group and normal peptide
 332 C group in case of shielding parameters
 333 c      write (iout,*) "rpp(1,1)",rpp(1,1)," pi",pi
 334 c      VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
 335 c      VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
 336 c      write (iout,*) "VSolvSphere",VSolvSphere,"VSolvSphere_div",
 337 c     &  VSolvSphere_div
 338 C long axis of side chain
 339 c      do i=1,ntyp
 340 c      long_r_sidechain(i)=vbldsc0(1,i)
 341 c      short_r_sidechain(i)=sigma0(i)
 342 c      enddo
 343 c      buff_shield=1.0d0
 344 c      endif
 345       if (me.eq.king .or. .not.out1file )
 346      &  write (iout,*) "DISTCHAINMAX",distchainmax
 347
 348       if(me.eq.king.or..not.out1file)
 349      & write (iout,'(2a)') diagmeth(kdiag),
 350      &  ' routine used to diagonalize matrices.'
 351       return
 352       end
 353 c--------------------------------------------------------------------------
 354       subroutine read_REMDpar
 355 C
 356 C Read REMD settings
 357 C
 358       implicit real*8 (a-h,o-z)
 359       include 'DIMENSIONS'
 360       include 'COMMON.IOUNITS'
 361       include 'COMMON.TIME1'
 362       include 'COMMON.MD'
 363 #ifndef LANG0
 364       include 'COMMON.LANGEVIN'
 365 #else
 366       include 'COMMON.LANGEVIN.lang0'
 367 #endif
 368       include 'COMMON.INTERACT'
 369       include 'COMMON.NAMES'
 370       include 'COMMON.GEO'
 371       include 'COMMON.REMD'
 372       include 'COMMON.CONTROL'
 373       include 'COMMON.SETUP'
 374       character*80 ucase
 375       character*320 controlcard
 376       character*3200 controlcard1
 377       integer iremd_m_total
 378
 379       if(me.eq.king.or..not.out1file)
 380      & write (iout,*) "REMD setup"
 381
 382       call card_concat(controlcard)
 383       call readi(controlcard,"NREP",nrep,3)
 384       call readi(controlcard,"NSTEX",nstex,1000)
 385       call reada(controlcard,"RETMIN",retmin,10.0d0)
 386       call reada(controlcard,"RETMAX",retmax,1000.0d0)
 387       mremdsync=(index(controlcard,'SYNC').gt.0)
 388       call readi(controlcard,"NSYN",i_sync_step,100)
 389       restart1file=(index(controlcard,'REST1FILE').gt.0)
 390       traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
 391       call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
 392       if(max_cache_traj_use.gt.max_cache_traj)
 393      &           max_cache_traj_use=max_cache_traj
 394       if(me.eq.king.or..not.out1file) then
 395 cd       if (traj1file) then
 396 crc caching is in testing - NTWX is not ignored
 397 cd        write (iout,*) "NTWX value is ignored"
 398 cd        write (iout,*) "  trajectory is stored to one file by master"
 399 cd        write (iout,*) "  before exchange at NSTEX intervals"
 400 cd       endif
 401        write (iout,*) "NREP= ",nrep
 402        write (iout,*) "NSTEX= ",nstex
 403        write (iout,*) "SYNC= ",mremdsync
 404        write (iout,*) "NSYN= ",i_sync_step
 405        write (iout,*) "TRAJCACHE= ",max_cache_traj_use
 406       endif
 407       remd_tlist=.false.
 408       if (index(controlcard,'TLIST').gt.0) then
 409          remd_tlist=.true.
 410          call card_concat(controlcard1)
 411          read(controlcard1,*) (remd_t(i),i=1,nrep)
 412          if(me.eq.king.or..not.out1file)
 413      &    write (iout,*)'tlist',(remd_t(i),i=1,nrep)
 414       endif
 415       remd_mlist=.false.
 416       if (index(controlcard,'MLIST').gt.0) then
 417          remd_mlist=.true.
 418          call card_concat(controlcard1)
 419          read(controlcard1,*) (remd_m(i),i=1,nrep)
 420          if(me.eq.king.or..not.out1file) then
 421           write (iout,*)'mlist',(remd_m(i),i=1,nrep)
 422           iremd_m_total=0
 423           do i=1,nrep
 424            iremd_m_total=iremd_m_total+remd_m(i)
 425           enddo
 426            write (iout,*) 'Total number of replicas ',iremd_m_total
 427           endif
 428          endif
 429         if (adaptive) then
 430           write (iout,*)
 431      &      "Adaptive (PMF-biased) umbrella sampling will be run"
 432           call PMFread
 433         endif
 434       if(me.eq.king.or..not.out1file)
 435      &   write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
 436       return
 437       end
 438 c--------------------------------------------------------------------------
 439       subroutine read_MDpar
 440 C
 441 C Read MD settings
 442 C
 443       implicit real*8 (a-h,o-z)
 444       include 'DIMENSIONS'
 445       include 'COMMON.IOUNITS'
 446       include 'COMMON.TIME1'
 447       include 'COMMON.MD'
 448 #ifndef LANG0
 449       include 'COMMON.LANGEVIN'
 450 #else
 451       include 'COMMON.LANGEVIN.lang0'
 452 #endif
 453       include 'COMMON.INTERACT'
 454       include 'COMMON.NAMES'
 455       include 'COMMON.GEO'
 456       include 'COMMON.SETUP'
 457       include 'COMMON.CONTROL'
 458       include 'COMMON.SPLITELE'
 459       include 'COMMON.FFIELD'
 460       character*80 ucase
 461       character*320 controlcard
 462
 463       call card_concat(controlcard)
 464       call readi(controlcard,"NSTEP",n_timestep,1000000)
 465       call readi(controlcard,"NTWE",ntwe,100)
 466       call readi(controlcard,"NTWX",ntwx,1000)
 467       call reada(controlcard,"DT",d_time,1.0d-1)
 468       call reada(controlcard,"DVMAX",dvmax,2.0d1)
 469       call reada(controlcard,"DAMAX",damax,1.0d1)
 470       call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
 471       call readi(controlcard,"LANG",lang,0)
 472       RESPA = index(controlcard,"RESPA") .gt. 0
 473       call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
 474       ntime_split0=ntime_split
 475       call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
 476       ntime_split0=ntime_split
 477 c      call reada(controlcard,"R_CUT",r_cut,2.0d0)
 478 c      call reada(controlcard,"LAMBDA",rlamb,0.3d0)
 479       rest = index(controlcard,"REST").gt.0
 480       tbf = index(controlcard,"TBF").gt.0
 481       usampl = index(controlcard,"USAMPL").gt.0
 482       scale_umb =  index(controlcard,"SCALE_UMB").gt.0
 483       adaptive = index(controlcard,"ADAPTIVE").gt.0
 484       mdpdb = index(controlcard,"MDPDB").gt.0
 485       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 486       call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
 487       call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
 488       call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
 489       if (count_reset_moment.eq.0) count_reset_moment=1000000000
 490       call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
 491       reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
 492       reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
 493       if (count_reset_vel.eq.0) count_reset_vel=1000000000
 494       large = index(controlcard,"LARGE").gt.0
 495       print_compon = index(controlcard,"PRINT_COMPON").gt.0
 496       rattle = index(controlcard,"RATTLE").gt.0
 497 c  if performing umbrella sampling, fragments constrained are read from the fragment file
 498       nset=0
 499       if(usampl) then
 500         write (iout,*) "Umbrella sampling will be run"
 501         if (scale_umb.and.adaptive) then
 502           write (iout,*) "ADAPTIVE and SCALE_UMB are mutually exclusive"
 503           write (iout,*) "Select one of those and re-run the job."
 504           stop
 505         endif
 506         if (scale_umb) write (iout,*)
 507      &"Umbrella-restraint force constants will be scaled by temperature"
 508         call read_fragments
 509       endif
 510
 511       if(me.eq.king.or..not.out1file) then
 512        write (iout,*)
 513        write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
 514        write (iout,*)
 515        write (iout,'(a)') "The units are:"
 516        write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
 517        write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
 518      &  " acceleration: angstrom/(48.9 fs)**2"
 519        write (iout,'(a)') "energy: kcal/mol, temperature: K"
 520        write (iout,*)
 521        write (iout,'(a60,i10)') "Number of time steps:",n_timestep
 522        write (iout,'(a60,f10.5,a)')
 523      &  "Initial time step of numerical integration:",d_time,
 524      &  " natural units"
 525        write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
 526        if (RESPA) then
 527         write (iout,'(2a,i4,a)')
 528      &    "A-MTS algorithm used; initial time step for fast-varying",
 529      &    " short-range forces split into",ntime_split," steps."
 530         write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
 531      &   r_cut," lambda",rlamb
 532        endif
 533        write (iout,'(2a,f10.5)')
 534      &  "Maximum acceleration threshold to reduce the time step",
 535      &  "/increase split number:",damax
 536        write (iout,'(2a,f10.5)')
 537      &  "Maximum predicted energy drift to reduce the timestep",
 538      &  "/increase split number:",edriftmax
 539        write (iout,'(a60,f10.5)')
 540      & "Maximum velocity threshold to reduce velocities:",dvmax
 541        write (iout,'(a60,i10)') "Frequency of property output:",ntwe
 542        write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
 543        if (rattle) write (iout,'(a60)')
 544      &  "Rattle algorithm used to constrain the virtual bonds"
 545       endif
 546       reset_fricmat=1000
 547       if (lang.gt.0) then
 548         call reada(controlcard,"ETAWAT",etawat,0.8904d0)
 549         call reada(controlcard,"RWAT",rwat,1.4d0)
 550         call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
 551         surfarea=index(controlcard,"SURFAREA").gt.0
 552         call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
 553         if(me.eq.king.or..not.out1file)then
 554          write (iout,'(/a,$)') "Langevin dynamics calculation"
 555          if (lang.eq.1) then
 556           write (iout,'(a/)')
 557      &      " with direct integration of Langevin equations"
 558          else if (lang.eq.2) then
 559           write (iout,'(a/)') " with TINKER stochasic MD integrator"
 560          else if (lang.eq.3) then
 561           write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
 562          else if (lang.eq.4) then
 563           write (iout,'(a/)') " in overdamped mode"
 564          else
 565           write (iout,'(//a,i5)')
 566      &      "=========== ERROR: Unknown Langevin dynamics mode:",lang
 567           stop
 568          endif
 569          write (iout,'(a60,f10.5)') "Temperature:",t_bath
 570          write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
 571          write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
 572          write (iout,'(a60,f10.5)')
 573      &   "Scaling factor of the friction forces:",scal_fric
 574          if (surfarea) write (iout,'(2a,i10,a)')
 575      &     "Friction coefficients will be scaled by solvent-accessible",
 576      &     " surface area every",reset_fricmat," steps."
 577         endif
 578 c Calculate friction coefficients and bounds of stochastic forces
 579         eta=6*pi*cPoise*etawat
 580         if(me.eq.king.or..not.out1file)
 581      &   write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
 582      &   ,eta
 583         gamp=scal_fric*(pstok+rwat)*eta
 584         stdfp=dsqrt(2*Rb*t_bath/d_time)
 585         do i=1,ntyp
 586           gamsc(i)=scal_fric*(restok(i)+rwat)*eta
 587           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
 588         enddo
 589         if(me.eq.king.or..not.out1file)then
 590          write (iout,'(/2a/)')
 591      &   "Radii of site types and friction coefficients and std's of",
 592      &   " stochastic forces of fully exposed sites"
 593          write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
 594          do i=1,ntyp
 595           write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
 596      &     gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
 597          enddo
 598         endif
 599       else if (tbf) then
 600         if(me.eq.king.or..not.out1file)then
 601          write (iout,'(a)') "Berendsen bath calculation"
 602          write (iout,'(a60,f10.5)') "Temperature:",t_bath
 603          write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
 604          if (reset_moment)
 605      &   write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
 606      &   count_reset_moment," steps"
 607          if (reset_vel)
 608      &    write (iout,'(a,i10,a)')
 609      &    "Velocities will be reset at random every",count_reset_vel,
 610      &   " steps"
 611         endif
 612       else
 613         if(me.eq.king.or..not.out1file)
 614      &   write (iout,'(a31)') "Microcanonical mode calculation"
 615       endif
 616       if(me.eq.king.or..not.out1file)then
 617        if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
 618        if (usampl) then
 619           write(iout,*) "MD running with constraints."
 620           write(iout,*) "Equilibration time ", eq_time, " mtus."
 621           write(iout,*) "Constraining ", nfrag," fragments."
 622           write(iout,*) "Length of each fragment, weight and q0:"
 623           do iset=1,nset
 624            write (iout,*) "Set of restraints #",iset
 625            do i=1,nfrag
 626               write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
 627      &           ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
 628            enddo
 629            write(iout,*) "constraints between ", npair, "fragments."
 630            write(iout,*) "constraint pairs, weights and q0:"
 631            do i=1,npair
 632             write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
 633      &             ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
 634            enddo
 635            write(iout,*) "angle constraints within ", nfrag_back,
 636      &      "backbone fragments."
 637            if (loc_qlike) then
 638            write(iout,*) "fragment, weights, q0:"
 639            do i=1,nfrag_back
 640             write(iout,'(2i5,3(f8.1,f8.2))') ifrag_back(1,i,iset),
 641      &         ifrag_back(2,i,iset),
 642      &         wfrag_back(1,i,iset),qin_back(1,i,iset),
 643      &         wfrag_back(2,i,iset),qin_back(2,i,iset),
 644      &         wfrag_back(3,i,iset),qin_back(3,i,iset)
 645            enddo
 646            else
 647            write(iout,*) "fragment, weights:"
 648            do i=1,nfrag_back
 649             write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
 650      &         ifrag_back(2,i,iset),wfrag_back(1,i,iset),
 651      &         wfrag_back(2,i,iset),wfrag_back(3,i,iset)
 652            enddo
 653            endif
 654           enddo
 655         iset=mod(kolor,nset)+1
 656        endif
 657       endif
 658       if(me.eq.king.or..not.out1file)
 659      & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
 660       return
 661       end
 662 c------------------------------------------------------------------------------
 663       subroutine molread
 664 C
 665 C Read molecular data.
 666 C
 667       implicit real*8 (a-h,o-z)
 668       include 'DIMENSIONS'
 669 #ifdef MPI
 670       include 'mpif.h'
 671       integer error_msg
 672 #endif
 673       include 'COMMON.IOUNITS'
 674       include 'COMMON.GEO'
 675       include 'COMMON.VAR'
 676       include 'COMMON.INTERACT'
 677       include 'COMMON.LOCAL'
 678       include 'COMMON.NAMES'
 679       include 'COMMON.CHAIN'
 680       include 'COMMON.FFIELD'
 681       include 'COMMON.SBRIDGE'
 682       include 'COMMON.HEADER'
 683       include 'COMMON.CONTROL'
 684       include 'COMMON.DBASE'
 685       include 'COMMON.THREAD'
 686       include 'COMMON.CONTACTS'
 687       include 'COMMON.TORCNSTR'
 688       include 'COMMON.TIME1'
 689       include 'COMMON.BOUNDS'
 690       include 'COMMON.MD'
 691       include 'COMMON.SETUP'
 692       include 'COMMON.SHIELD'
 693       character*4 sequence(maxres)
 694       integer rescode
 695       double precision x(maxvar)
 696       character*256 pdbfile
 697       character*400 weightcard
 698       character*80 weightcard_t,ucase
 699       dimension itype_pdb(maxres)
 700       common /pizda/ itype_pdb
 701       logical seq_comp,fail
 702       double precision energia(0:n_ene)
 703       double precision secprob(3,maxdih_constr)
 704       integer ilen
 705       external ilen
 706 C
 707 C Read PDB structure if applicable
 708 C
 709       if (indpdb.gt.0 .or. pdbref) then
 710         read(inp,'(a)') pdbfile
 711         if(me.eq.king.or..not.out1file)
 712      &   write (iout,'(2a)') 'PDB data will be read from file ',
 713      &   pdbfile(:ilen(pdbfile))
 714         open(ipdbin,file=pdbfile,status='old',err=33)
 715         goto 34
 716   33    write (iout,'(a)') 'Error opening PDB file.'
 717         stop
 718   34    continue
 719         call readpdb
 720         if(me.eq.king.or..not.out1file)
 721      &   write (iout,'(a,i3,a,i3)')'nsup=',nsup,
 722      &   ' nstart_sup=',nstart_sup
 723         do i=1,nres
 724           itype_pdb(i)=itype(i)
 725         enddo
 726         close (ipdbin)
 727         nnt=nstart_sup
 728         nct=nstart_sup+nsup-1
 729         call contact(.false.,ncont_ref,icont_ref,co)
 730
 731         if (sideadd) then
 732          if(me.eq.king.or..not.out1file)
 733      &    write(iout,*)'Adding sidechains'
 734          maxsi=1000
 735          do i=2,nres-1
 736           iti=itype(i)
 737           if (iti.ne.10 .and. itype(i).ne.ntyp1) then
 738             nsi=0
 739             fail=.true.
 740             do while (fail.and.nsi.le.maxsi)
 741               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
 742               nsi=nsi+1
 743             enddo
 744             if(fail) write(iout,*)'Adding sidechain failed for res ',
 745      &              i,' after ',nsi,' trials'
 746           endif
 747          enddo
 748         endif
 749       endif
 750       if (indpdb.eq.0) then
 751 C Read sequence if not taken from the pdb file.
 752         read (inp,*) nres
 753 c        print *,'nres=',nres
 754         if (iscode.gt.0) then
 755           read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 756         else
 757           read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 758         endif
 759 C Convert sequence to numeric code
 760         do i=1,nres
 761           itype(i)=rescode(i,sequence(i),iscode)
 762         enddo
 763 C Assign initial virtual bond lengths
 764         do i=2,nres
 765           vbld(i)=vbl
 766           vbld_inv(i)=vblinv
 767         enddo
 768         do i=2,nres-1
 769           vbld(i+nres)=dsc(iabs(itype(i)))
 770           vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
 771 c          write (iout,*) "i",i," itype",itype(i),
 772 c     &      " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
 773         enddo
 774       endif
 775 c      print *,nres
 776 c      print '(20i4)',(itype(i),i=1,nres)
 777       do i=1,nres
 778 #ifdef PROCOR
 779         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
 780 #else
 781         if (itype(i).eq.ntyp1) then
 782 #endif
 783           itel(i)=0
 784 #ifdef PROCOR
 785         else if (iabs(itype(i+1)).ne.20) then
 786 #else
 787         else if (iabs(itype(i)).ne.20) then
 788 #endif
 789           itel(i)=1
 790         else
 791           itel(i)=2
 792         endif
 793       enddo
 794       if(me.eq.king.or..not.out1file)then
 795        write (iout,*) "ITEL"
 796        do i=1,nres-1
 797          write (iout,*) i,itype(i),itel(i)
 798        enddo
 799        print *,'Call Read_Bridge.'
 800       endif
 801       nnt=1
 802       nct=nres
 803 cd      print *,'NNT=',NNT,' NCT=',NCT
 804       call seq2chains(nres,itype,nchain,chain_length,chain_border,
 805      &  ireschain)
 806       write(iout,*) "nres",nres," nchain",nchain
 807       do i=1,nchain
 808         write(iout,*)"chain",i,chain_length(i),chain_border(1,i),
 809      &    chain_border(2,i)
 810       enddo
 811       call chain_symmetry(nchain,nres,itype,chain_border,
 812      &    chain_length,npermchain,tabpermchain)
 813       do i=1,nres
 814         write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain),
 815      &    ii=1,npermchain)
 816       enddo
 817       write(iout,*) "residue permutations"
 818       do i=1,nres
 819         write(iout,*) i,(iperm(i,ii),ii=1,npermchain)
 820       enddo
 821       if (itype(1).eq.ntyp1) nnt=2
 822       if (itype(nres).eq.ntyp1) nct=nct-1
 823       if (pdbref) then
 824         if(me.eq.king.or..not.out1file)
 825      &   write (iout,'(a,i3)') 'nsup=',nsup
 826         nstart_seq=nnt
 827         if (nsup.le.(nct-nnt+1)) then
 828           do i=0,nct-nnt+1-nsup
 829             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 830               nstart_seq=nnt+i
 831               goto 111
 832             endif
 833           enddo
 834           write (iout,'(a)')
 835      &            'Error - sequences to be superposed do not match.'
 836           stop
 837         else
 838           do i=0,nsup-(nct-nnt+1)
 839             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 840      &      then
 841               nstart_sup=nstart_sup+i
 842               nsup=nct-nnt+1
 843               goto 111
 844             endif
 845           enddo
 846           write (iout,'(a)')
 847      &            'Error - sequences to be superposed do not match.'
 848         endif
 849   111   continue
 850         if (nsup.eq.0) nsup=nct-nnt
 851         if (nstart_sup.eq.0) nstart_sup=nnt
 852         if (nstart_seq.eq.0) nstart_seq=nnt
 853         if(me.eq.king.or..not.out1file)
 854      &   write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
 855      &                 ' nstart_seq=',nstart_seq
 856       endif
 857       call read_bridge
 858 C 8/13/98 Set limits to generating the dihedral angles
 859       do i=1,nres
 860         phibound(1,i)=-pi
 861         phibound(2,i)=pi
 862       enddo
 863       read (inp,*) ndih_constr
 864       write (iout,*) "ndih_constr",ndih_constr
 865       if (ndih_constr.gt.0) then
 866         raw_psipred=.false.
 867 C        read (inp,*) ftors
 868         read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
 869      &  i=1,ndih_constr)
 870 #ifdef MPI
 871         if(me.eq.king.or..not.out1file)then
 872 #endif
 873          write (iout,*)
 874      &   'There are',ndih_constr,' restraints on gamma angles.'
 875          do i=1,ndih_constr
 876           write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
 877      &    ftors(i)
 878          enddo
 879 #ifdef MPI
 880         endif
 881 #endif
 882         do i=1,ndih_constr
 883           phi0(i)=deg2rad*phi0(i)
 884           drange(i)=deg2rad*drange(i)
 885         enddo
 886 C        if(me.eq.king.or..not.out1file)
 887 C     &   write (iout,*) 'FTORS',ftors
 888         do i=1,ndih_constr
 889           ii = idih_constr(i)
 890           phibound(1,ii) = phi0(i)-drange(i)
 891           phibound(2,ii) = phi0(i)+drange(i)
 892         enddo
 893 #ifdef MPI
 894         if (me.eq.king .or. .not.out1file) then
 895 #endif
 896          write (iout,'(a)') 'Boundaries in gamma angle sampling:'
 897          do i=1,nres
 898            write (iout,'(a3,i5,2f10.1)')
 899      &    restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
 900          enddo
 901 #ifdef MP
 902         endif
 903 #endif
 904       else if (ndih_constr.lt.0) then
 905         raw_psipred=.true.
 906         call card_concat(weightcard)
 907         call reada(weightcard,"PHIHEL",phihel,50.0D0)
 908         call reada(weightcard,"PHIBET",phibet,180.0D0)
 909         call reada(weightcard,"SIGMAHEL",sigmahel,30.0d0)
 910         call reada(weightcard,"SIGMABET",sigmabet,40.0d0)
 911         call reada(weightcard,"WDIHC",wdihc,0.591D0)
 912         write (iout,*) "Weight of dihedral angle restraints",wdihc
 913         read(inp,'(9x,3f7.3)')
 914      &     (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
 915         write (iout,*) "The secprob array"
 916         do i=nnt,nct
 917           write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
 918         enddo
 919         ndih_constr=0
 920         do i=nnt+3,nct
 921           if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
 922      &    .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
 923             ndih_constr=ndih_constr+1
 924             idih_constr(ndih_constr)=i
 925             sumv=0.0d0
 926             do j=1,3
 927               vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
 928               sumv=sumv+vpsipred(j,ndih_constr)
 929             enddo
 930             do j=1,3
 931               vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
 932             enddo
 933             phibound(1,ndih_constr)=phihel*deg2rad
 934             phibound(2,ndih_constr)=phibet*deg2rad
 935             sdihed(1,ndih_constr)=sigmahel*deg2rad
 936             sdihed(2,ndih_constr)=sigmabet*deg2rad
 937           endif
 938         enddo
 939 #ifdef MPI
 940         if(me.eq.king.or..not.out1file)then
 941 #endif
 942          write (iout,*)
 943      &   'There are',ndih_constr,
 944      &   ' bimodal restraints on gamma angles.'
 945          do i=1,ndih_constr
 946            write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,
 947      &      restyp(itype(idih_constr(i)-2)),idih_constr(i)-2,
 948      &      restyp(itype(idih_constr(i)-1)),idih_constr(i)-1,
 949      &      phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,
 950      &      phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,
 951      &      (vpsipred(j,i),j=1,3)
 952          enddo
 953 #ifdef MPI
 954         endif
 955 #endif
 956 c Raw psipred input
 957       endif
 958 C first setting the theta boundaries to 0 to pi
 959 C this mean that there is no energy penalty for any angle occuring this can be applied
 960 C for generate random conformation but is not implemented in this way
 961 C      do i=1,nres
 962 C        thetabound(1,i)=0
 963 C        thetabound(2,i)=pi
 964 C      enddo
 965 C begin reading theta constrains this is quartic constrains allowing to
 966 C have smooth second derivative
 967       if (with_theta_constr) then
 968 C with_theta_constr is keyword allowing for occurance of theta constrains
 969       read (inp,*) ntheta_constr
 970 C ntheta_constr is the number of theta constrains
 971       if (ntheta_constr.gt.0) then
 972 C        read (inp,*) ftors
 973         read (inp,*) (itheta_constr(i),theta_constr0(i),
 974      &  theta_drange(i),for_thet_constr(i),
 975      &  i=1,ntheta_constr)
 976 C the above code reads from 1 to ntheta_constr
 977 C itheta_constr(i) residue i for which is theta_constr
 978 C theta_constr0 the global minimum value
 979 C theta_drange is range for which there is no energy penalty
 980 C for_thet_constr is the force constant for quartic energy penalty
 981 C E=k*x**4
 982         if(me.eq.king.or..not.out1file)then
 983          write (iout,*)
 984      &   'There are',ntheta_constr,' constraints on phi angles.'
 985          do i=1,ntheta_constr
 986           write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
 987      &    theta_drange(i),
 988      &    for_thet_constr(i)
 989          enddo
 990         endif
 991         do i=1,ntheta_constr
 992           theta_constr0(i)=deg2rad*theta_constr0(i)
 993           theta_drange(i)=deg2rad*theta_drange(i)
 994         enddo
 995 C        if(me.eq.king.or..not.out1file)
 996 C     &   write (iout,*) 'FTORS',ftors
 997 C        do i=1,ntheta_constr
 998 C          ii = itheta_constr(i)
 999 C          thetabound(1,ii) = phi0(i)-drange(i)
1000 C          thetabound(2,ii) = phi0(i)+drange(i)
1001 C        enddo
1002       endif ! ntheta_constr.gt.0
1003       endif! with_theta_constr
1004 C
1005 C      with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
1006 C      write (iout,*) "with_dihed_constr ",with_dihed_constr
1007 c--- Zscore rms -------
1008       if (nz_start.eq.0) nz_start=nnt
1009       if (nz_end.eq.0 .and. nsup.gt.0) then
1010         nz_end=nnt+nsup-1
1011       else if (nz_end.eq.0) then
1012         nz_end=nct
1013       endif
1014       if(me.eq.king.or..not.out1file)then
1015        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1016        write (iout,*) 'IZ_SC=',iz_sc
1017       endif
1018 c----------------------
1019       call init_int_table
1020       if (refstr) then
1021         if (.not.pdbref) then
1022           call read_angles(inp,*38)
1023           goto 39
1024    38     write (iout,'(a)') 'Error reading reference structure.'
1025 #ifdef MPI
1026           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1027           stop 'Error reading reference structure'
1028 #endif
1029    39     call chainbuild_extconf
1030           call setup_var
1031 czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
1032           nstart_sup=nnt
1033           nstart_seq=nnt
1034           nsup=nct-nnt+1
1035           do i=1,2*nres
1036             do j=1,3
1037               cref(j,i)=c(j,i)
1038             enddo
1039           enddo
1040           call contact(.true.,ncont_ref,icont_ref,co)
1041         endif
1042       endif
1043 c        print *, "A TU"
1044       write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1045       call flush(iout)
1046       if (constr_dist.gt.0) call read_dist_constr
1047       write (iout,*) "After read_dist_constr nhpb",nhpb
1048       if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1049       call hpb_partition
1050       call NMRpeak_partition
1051       if(me.eq.king.or..not.out1file)write (iout,*) 'Contact order:',co
1052       if (pdbref) then
1053         if(me.eq.king.or..not.out1file)
1054      &   write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1055         do i=1,ncont_ref
1056           do j=1,2
1057             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1058           enddo
1059           if(me.eq.king.or..not.out1file)
1060      &     write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1061      &     icont_ref(1,i),' ',
1062      &     restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1063         enddo
1064       endif
1065       write (iout,*) "calling read_saxs_consrtr",nsaxs
1066       if (nsaxs.gt.0) call read_saxs_constr
1067
1068 C      endif
1069       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1070      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1071      &    modecalc.ne.10) then
1072 C If input structure hasn't been supplied from the PDB file read or generate
1073 C initial geometry.
1074         if (iranconf.eq.0 .and. .not. extconf) then
1075           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1076      &     write (iout,'(a)') 'Initial geometry will be read in.'
1077           if (read_cart) then
1078             read(inp,'(8f10.5)',end=36,err=36)
1079      &       ((c(l,k),l=1,3),k=1,nres),
1080      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
1081             write (iout,*) "Exit READ_CART"
1082 c            write (iout,'(8f10.5)')
1083 c     &       ((c(l,k),l=1,3),k=1,nres),
1084 c     &       ((c(l,k+nres),l=1,3),k=nnt,nct)
1085             call cartprint
1086             do i=1,nres-1
1087               do j=1,3
1088                 dc(j,i)=c(j,i+1)-c(j,i)
1089 c                dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1090               enddo
1091             enddo
1092             do i=nnt,nct
1093               if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1094                 do j=1,3
1095                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
1096 c                  dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1097                 enddo
1098               else
1099                 do j=1,3
1100                   dc(j,i+nres)=0.0d0
1101 c                  dc_norm(j,i+nres)=0.0d0
1102                 enddo
1103               endif
1104             enddo
1105             call int_from_cart1(.true.)
1106             write (iout,*) "Finish INT_TO_CART"
1107 c      write (iout,*) "DC"
1108 c      do i=1,nres
1109 c        write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1110 c     &    (dc(j,i+nres),j=1,3)
1111 c      enddo
1112 c            call cartprint
1113 c            call setup_var
1114 c            return
1115           else
1116             call read_angles(inp,*36)
1117             call chainbuild_extconf
1118           endif
1119           goto 37
1120    36     write (iout,'(a)') 'Error reading angle file.'
1121 #ifdef MPI
1122           call mpi_finalize( MPI_COMM_WORLD,IERR )
1123 #endif
1124           stop 'Error reading angle file.'
1125    37     continue
1126         else if (extconf) then
1127          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1128      &    write (iout,'(a)') 'Extended chain initial geometry.'
1129          do i=3,nres
1130           theta(i)=90d0*deg2rad
1131          enddo
1132          do i=4,nres
1133           phi(i)=180d0*deg2rad
1134          enddo
1135          do i=2,nres-1
1136           alph(i)=110d0*deg2rad
1137          enddo
1138          do i=2,nres-1
1139           omeg(i)=-120d0*deg2rad
1140           if (itype(i).le.0) omeg(i)=-omeg(i)
1141          enddo
1142         call chainbuild_extconf
1143         else
1144           if(me.eq.king.or..not.out1file)
1145      &     write (iout,'(a)') 'Random-generated initial geometry.'
1146
1147
1148 #ifdef MPI
1149           if (me.eq.king  .or. fg_rank.eq.0 .and. (
1150      &           modecalc.eq.12 .or. modecalc.eq.14) ) then
1151 #endif
1152             do itrial=1,100
1153               itmp=1
1154               call gen_rand_conf(itmp,*30)
1155               goto 40
1156    30         write (iout,*) 'Failed to generate random conformation',
1157      &          ', itrial=',itrial
1158               write (*,*) 'Processor:',me,
1159      &          ' Failed to generate random conformation',
1160      &          ' itrial=',itrial
1161               call intout
1162
1163 #ifdef AIX
1164               call flush_(iout)
1165 #else
1166               call flush(iout)
1167 #endif
1168             enddo
1169             write (iout,'(a,i3,a)') 'Processor:',me,
1170      &        ' error in generating random conformation.'
1171             write (*,'(a,i3,a)') 'Processor:',me,
1172      &        ' error in generating random conformation.'
1173             call flush(iout)
1174 #ifdef MPI
1175             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1176    40       continue
1177           endif
1178 #else
1179           write (*,'(a)')
1180      &      ' error in generating random conformation.'
1181           stop
1182    40     continue
1183 #endif
1184         endif
1185       elseif (modecalc.eq.4) then
1186         read (inp,'(a)') intinname
1187         open (intin,file=intinname,status='old',err=333)
1188         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1189      &  write (iout,'(a)') 'intinname',intinname
1190         write (*,'(a)') 'Processor',myrank,' intinname',intinname
1191         goto 334
1192   333   write (iout,'(2a)') 'Error opening angle file ',intinname
1193 #ifdef MPI
1194         call MPI_Finalize(MPI_COMM_WORLD,IERR)
1195 #endif
1196         stop 'Error opening angle file.'
1197   334   continue
1198
1199       endif
1200 C Generate distance constraints, if the PDB structure is to be regularized.
1201       if (nthread.gt.0) then
1202         call read_threadbase
1203       endif
1204       call setup_var
1205       if (me.eq.king .or. .not. out1file)
1206      & call intout
1207       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1208         write (iout,'(/a,i3,a)')
1209      &  'The chain contains',ns,' disulfide-bridging cysteines.'
1210         write (iout,'(20i4)') (iss(i),i=1,ns)
1211        if (dyn_ss) then
1212           write(iout,*)"Running with dynamic disulfide-bond formation"
1213        else
1214         write (iout,'(/a/)') 'Pre-formed links are:'
1215         do i=1,nss
1216           i1=ihpb(i)-nres
1217           i2=jhpb(i)-nres
1218           it1=itype(i1)
1219           it2=itype(i2)
1220           write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1221      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1222      &    ebr,forcon(i)
1223         enddo
1224         write (iout,'(a)')
1225        endif
1226       endif
1227       if (ns.gt.0.and.dyn_ss) then
1228           do i=nss+1,nhpb
1229             ihpb(i-nss)=ihpb(i)
1230             jhpb(i-nss)=jhpb(i)
1231             forcon(i-nss)=forcon(i)
1232             dhpb(i-nss)=dhpb(i)
1233           enddo
1234           nhpb=nhpb-nss
1235           nss=0
1236           call hpb_partition
1237           do i=1,ns
1238             dyn_ss_mask(iss(i))=.true.
1239           enddo
1240       endif
1241 c --- test
1242 c      call cartprint
1243 c      write (iout,*) "DC"
1244 c      do i=1,nres
1245 c        write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1246 c     &    (dc(j,i+nres),j=1,3)
1247 c      enddo
1248       if (i2ndstr.gt.0) call secstrp2dihc
1249 c      call geom_to_var(nvar,x)
1250 c      call etotal(energia(0))
1251 c      call enerprint(energia(0))
1252 c      call briefout(0,etot)
1253 c      stop
1254 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1255 cd    write (iout,'(a)') 'Variable list:'
1256 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1257 #ifdef MPI
1258       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1259      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1260      &  'Processor',myrank,': end reading molecular data.'
1261 #endif
1262 c      print *,"A TU?"
1263       return
1264       end
1265 c--------------------------------------------------------------------------
1266       logical function seq_comp(itypea,itypeb,length)
1267       implicit none
1268       integer length,itypea(length),itypeb(length)
1269       integer i
1270       do i=1,length
1271         if (itypea(i).ne.itypeb(i)) then
1272           seq_comp=.false.
1273           return
1274         endif
1275       enddo
1276       seq_comp=.true.
1277       return
1278       end
1279 c-----------------------------------------------------------------------------
1280       subroutine read_bridge
1281 C Read information about disulfide bridges.
1282       implicit real*8 (a-h,o-z)
1283       include 'DIMENSIONS'
1284 #ifdef MPI
1285       include 'mpif.h'
1286 #endif
1287       include 'COMMON.IOUNITS'
1288       include 'COMMON.GEO'
1289       include 'COMMON.VAR'
1290       include 'COMMON.INTERACT'
1291       include 'COMMON.LOCAL'
1292       include 'COMMON.NAMES'
1293       include 'COMMON.CHAIN'
1294       include 'COMMON.FFIELD'
1295       include 'COMMON.SBRIDGE'
1296       include 'COMMON.HEADER'
1297       include 'COMMON.CONTROL'
1298       include 'COMMON.DBASE'
1299       include 'COMMON.THREAD'
1300       include 'COMMON.TIME1'
1301       include 'COMMON.SETUP'
1302 C Read bridging residues.
1303       read (inp,*) ns,(iss(i),i=1,ns)
1304       print *,'ns=',ns
1305       if(me.eq.king.or..not.out1file)
1306      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1307 C Check whether the specified bridging residues are cystines.
1308       do i=1,ns
1309         if (itype(iss(i)).ne.1) then
1310           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1311      &   'Do you REALLY think that the residue ',
1312      &    restyp(itype(iss(i))),i,
1313      &   ' can form a disulfide bridge?!!!'
1314           write (*,'(2a,i3,a)')
1315      &   'Do you REALLY think that the residue ',
1316      &    restyp(itype(iss(i))),i,
1317      &   ' can form a disulfide bridge?!!!'
1318 #ifdef MPI
1319          call MPI_Finalize(MPI_COMM_WORLD,ierror)
1320          stop
1321 #endif
1322         endif
1323       enddo
1324 C Read preformed bridges.
1325       if (ns.gt.0) then
1326       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1327       if(fg_rank.eq.0)
1328      & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1329       if (nss.gt.0) then
1330         nhpb=nss
1331 C Check if the residues involved in bridges are in the specified list of
1332 C bridging residues.
1333         do i=1,nss
1334           do j=1,i-1
1335             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1336      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1337               write (iout,'(a,i3,a)') 'Disulfide pair',i,
1338      &      ' contains residues present in other pairs.'
1339               write (*,'(a,i3,a)') 'Disulfide pair',i,
1340      &      ' contains residues present in other pairs.'
1341 #ifdef MPI
1342               call MPI_Finalize(MPI_COMM_WORLD,ierror)
1343               stop
1344 #endif
1345             endif
1346           enddo
1347           do j=1,ns
1348             if (ihpb(i).eq.iss(j)) goto 10
1349           enddo
1350           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1351    10     continue
1352           do j=1,ns
1353             if (jhpb(i).eq.iss(j)) goto 20
1354           enddo
1355           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1356    20     continue
1357           dhpb(i)=dbr
1358           forcon(i)=fbr
1359         enddo
1360         do i=1,nss
1361           ihpb(i)=ihpb(i)+nres
1362           jhpb(i)=jhpb(i)+nres
1363         enddo
1364       endif
1365       endif
1366       return
1367       end
1368 c----------------------------------------------------------------------------
1369       subroutine read_x(kanal,*)
1370       implicit real*8 (a-h,o-z)
1371       include 'DIMENSIONS'
1372       include 'COMMON.GEO'
1373       include 'COMMON.VAR'
1374       include 'COMMON.CHAIN'
1375       include 'COMMON.IOUNITS'
1376       include 'COMMON.CONTROL'
1377       include 'COMMON.LOCAL'
1378       include 'COMMON.INTERACT'
1379 c Read coordinates from input
1380 c
1381       read(kanal,'(8f10.5)',end=10,err=10)
1382      &  ((c(l,k),l=1,3),k=1,nres),
1383      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
1384       do j=1,3
1385         c(j,nres+1)=c(j,1)
1386         c(j,2*nres)=c(j,nres)
1387       enddo
1388       call int_from_cart1(.false.)
1389       do i=1,nres-1
1390         do j=1,3
1391           dc(j,i)=c(j,i+1)-c(j,i)
1392           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1393         enddo
1394       enddo
1395       do i=nnt,nct
1396         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1397           do j=1,3
1398             dc(j,i+nres)=c(j,i+nres)-c(j,i)
1399             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1400           enddo
1401         endif
1402       enddo
1403
1404       return
1405    10 return1
1406       end
1407 c----------------------------------------------------------------------------
1408       subroutine read_threadbase
1409       implicit real*8 (a-h,o-z)
1410       include 'DIMENSIONS'
1411       include 'COMMON.IOUNITS'
1412       include 'COMMON.GEO'
1413       include 'COMMON.VAR'
1414       include 'COMMON.INTERACT'
1415       include 'COMMON.LOCAL'
1416       include 'COMMON.NAMES'
1417       include 'COMMON.CHAIN'
1418       include 'COMMON.FFIELD'
1419       include 'COMMON.SBRIDGE'
1420       include 'COMMON.HEADER'
1421       include 'COMMON.CONTROL'
1422       include 'COMMON.DBASE'
1423       include 'COMMON.THREAD'
1424       include 'COMMON.TIME1'
1425 C Read pattern database for threading.
1426       read (icbase,*) nseq
1427       do i=1,nseq
1428         read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1429      &   nres_base(2,i),nres_base(3,i)
1430         read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1431      &   nres_base(1,i))
1432 c       write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1433 c    &   nres_base(2,i),nres_base(3,i)
1434 c       write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1435 c    &   nres_base(1,i))
1436       enddo
1437       close (icbase)
1438       if (weidis.eq.0.0D0) weidis=0.1D0
1439       do i=nnt,nct
1440         do j=i+2,nct
1441           nhpb=nhpb+1
1442           ihpb(nhpb)=i
1443           jhpb(nhpb)=j
1444           forcon(nhpb)=weidis
1445         enddo
1446       enddo
1447       read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1448       write (iout,'(a,i5)') 'nexcl: ',nexcl
1449       write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1450       return
1451       end
1452 c------------------------------------------------------------------------------
1453       subroutine setup_var
1454       implicit real*8 (a-h,o-z)
1455       include 'DIMENSIONS'
1456       include 'COMMON.IOUNITS'
1457       include 'COMMON.GEO'
1458       include 'COMMON.VAR'
1459       include 'COMMON.INTERACT'
1460       include 'COMMON.LOCAL'
1461       include 'COMMON.NAMES'
1462       include 'COMMON.CHAIN'
1463       include 'COMMON.FFIELD'
1464       include 'COMMON.SBRIDGE'
1465       include 'COMMON.HEADER'
1466       include 'COMMON.CONTROL'
1467       include 'COMMON.DBASE'
1468       include 'COMMON.THREAD'
1469       include 'COMMON.TIME1'
1470 C Set up variable list.
1471       ntheta=nres-2
1472       nphi=nres-3
1473       nvar=ntheta+nphi
1474       nside=0
1475       write (iout,*) "SETUP_VAR ialph"
1476       do i=2,nres-1
1477         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1478           nside=nside+1
1479           ialph(i,1)=nvar+nside
1480           ialph(nside,2)=i
1481         endif
1482       enddo
1483       if (indphi.gt.0) then
1484         nvar=nphi
1485       else if (indback.gt.0) then
1486         nvar=nphi+ntheta
1487       else
1488         nvar=nvar+2*nside
1489       endif
1490       write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1491       return
1492       end
1493 c----------------------------------------------------------------------------
1494       subroutine gen_dist_constr
1495 C Generate CA distance constraints.
1496       implicit real*8 (a-h,o-z)
1497       include 'DIMENSIONS'
1498       include 'COMMON.IOUNITS'
1499       include 'COMMON.GEO'
1500       include 'COMMON.VAR'
1501       include 'COMMON.INTERACT'
1502       include 'COMMON.LOCAL'
1503       include 'COMMON.NAMES'
1504       include 'COMMON.CHAIN'
1505       include 'COMMON.FFIELD'
1506       include 'COMMON.SBRIDGE'
1507       include 'COMMON.HEADER'
1508       include 'COMMON.CONTROL'
1509       include 'COMMON.DBASE'
1510       include 'COMMON.THREAD'
1511       include 'COMMON.TIME1'
1512       dimension itype_pdb(maxres)
1513       common /pizda/ itype_pdb
1514       character*2 iden
1515 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1516 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1517 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1518 cd     & ' nsup',nsup
1519       do i=nstart_sup,nstart_sup+nsup-1
1520 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1521 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1522         do j=i+2,nstart_sup+nsup-1
1523           nhpb=nhpb+1
1524           ihpb(nhpb)=i+nstart_seq-nstart_sup
1525           jhpb(nhpb)=j+nstart_seq-nstart_sup
1526           forcon(nhpb)=weidis
1527           dhpb(nhpb)=dist(i,j)
1528         enddo
1529       enddo
1530 cd      write (iout,'(a)') 'Distance constraints:'
1531 cd      do i=nss+1,nhpb
1532 cd        ii=ihpb(i)
1533 cd        jj=jhpb(i)
1534 cd        iden='CA'
1535 cd        if (ii.gt.nres) then
1536 cd          iden='SC'
1537 cd          ii=ii-nres
1538 cd          jj=jj-nres
1539 cd        endif
1540 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1541 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1542 cd     &  dhpb(i),forcon(i)
1543 cd      enddo
1544       return
1545       end
1546 c----------------------------------------------------------------------------
1547       subroutine map_read
1548       implicit real*8 (a-h,o-z)
1549       include 'DIMENSIONS'
1550       include 'COMMON.MAP'
1551       include 'COMMON.IOUNITS'
1552       character*3 angid(4) /'THE','PHI','ALP','OME'/
1553       character*80 mapcard,ucase
1554       do imap=1,nmap
1555         read (inp,'(a)') mapcard
1556         mapcard=ucase(mapcard)
1557         if (index(mapcard,'PHI').gt.0) then
1558           kang(imap)=1
1559         else if (index(mapcard,'THE').gt.0) then
1560           kang(imap)=2
1561         else if (index(mapcard,'ALP').gt.0) then
1562           kang(imap)=3
1563         else if (index(mapcard,'OME').gt.0) then
1564           kang(imap)=4
1565         else
1566           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1567           stop 'Error - illegal variable spec in MAP card.'
1568         endif
1569         call readi (mapcard,'RES1',res1(imap),0)
1570         call readi (mapcard,'RES2',res2(imap),0)
1571         if (res1(imap).eq.0) then
1572           res1(imap)=res2(imap)
1573         else if (res2(imap).eq.0) then
1574           res2(imap)=res1(imap)
1575         endif
1576         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1577           write (iout,'(a)')
1578      &    'Error - illegal definition of variable group in MAP.'
1579           stop 'Error - illegal definition of variable group in MAP.'
1580         endif
1581         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1582         call reada(mapcard,'TO',ang_to(imap),0.0D0)
1583         call readi(mapcard,'NSTEP',nstep(imap),0)
1584         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1585           write (iout,'(a)')
1586      &     'Illegal boundary and/or step size specification in MAP.'
1587           stop 'Illegal boundary and/or step size specification in MAP.'
1588         endif
1589       enddo ! imap
1590       return
1591       end
1592 c----------------------------------------------------------------------------
1593       subroutine csaread
1594       implicit real*8 (a-h,o-z)
1595       include 'DIMENSIONS'
1596       include 'COMMON.IOUNITS'
1597       include 'COMMON.GEO'
1598       include 'COMMON.CSA'
1599       include 'COMMON.BANK'
1600       include 'COMMON.CONTROL'
1601       character*80 ucase
1602       character*620 mcmcard
1603       call card_concat(mcmcard)
1604
1605       call readi(mcmcard,'NCONF',nconf,50)
1606       call readi(mcmcard,'NADD',nadd,0)
1607       call readi(mcmcard,'JSTART',jstart,1)
1608       call readi(mcmcard,'JEND',jend,1)
1609       call readi(mcmcard,'NSTMAX',nstmax,500000)
1610       call readi(mcmcard,'N0',n0,1)
1611       call readi(mcmcard,'N1',n1,6)
1612       call readi(mcmcard,'N2',n2,4)
1613       call readi(mcmcard,'N3',n3,0)
1614       call readi(mcmcard,'N4',n4,0)
1615       call readi(mcmcard,'N5',n5,0)
1616       call readi(mcmcard,'N6',n6,10)
1617       call readi(mcmcard,'N7',n7,0)
1618       call readi(mcmcard,'N8',n8,0)
1619       call readi(mcmcard,'N9',n9,0)
1620       call readi(mcmcard,'N14',n14,0)
1621       call readi(mcmcard,'N15',n15,0)
1622       call readi(mcmcard,'N16',n16,0)
1623       call readi(mcmcard,'N17',n17,0)
1624       call readi(mcmcard,'N18',n18,0)
1625
1626       vdisulf=(index(mcmcard,'DYNSS').gt.0)
1627
1628       call readi(mcmcard,'NDIFF',ndiff,2)
1629       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1630       call readi(mcmcard,'IS1',is1,1)
1631       call readi(mcmcard,'IS2',is2,8)
1632       call readi(mcmcard,'NRAN0',nran0,4)
1633       call readi(mcmcard,'NRAN1',nran1,2)
1634       call readi(mcmcard,'IRR',irr,1)
1635       call readi(mcmcard,'NSEED',nseed,20)
1636       call readi(mcmcard,'NTOTAL',ntotal,10000)
1637       call reada(mcmcard,'CUT1',cut1,2.0d0)
1638       call reada(mcmcard,'CUT2',cut2,5.0d0)
1639       call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1640       call readi(mcmcard,'ICMAX',icmax,3)
1641       call readi(mcmcard,'IRESTART',irestart,0)
1642 c!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
1643       ntbankm=0
1644 c!bankt
1645       call reada(mcmcard,'DELE',dele,20.0d0)
1646       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1647       call readi(mcmcard,'IREF',iref,0)
1648       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1649       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1650       call readi(mcmcard,'NCONF_IN',nconf_in,0)
1651       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1652       write (iout,*) "NCONF_IN",nconf_in
1653       return
1654       end
1655 c----------------------------------------------------------------------------
1656 cfmc      subroutine mcmfread
1657 cfmc      implicit real*8 (a-h,o-z)
1658 cfmc      include 'DIMENSIONS'
1659 cfmc      include 'COMMON.MCMF'
1660 cfmc      include 'COMMON.IOUNITS'
1661 cfmc      include 'COMMON.GEO'
1662 cfmc      character*80 ucase
1663 cfmc      character*620 mcmcard
1664 cfmc      call card_concat(mcmcard)
1665 cfmc
1666 cfmc      call readi(mcmcard,'MAXRANT',maxrant,1000)
1667 cfmc      write(iout,*)'MAXRANT=',maxrant
1668 cfmc      call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1669 cfmc      write(iout,*)'MAXFAM=',maxfam
1670 cfmc      call readi(mcmcard,'NNET1',nnet1,5)
1671 cfmc      write(iout,*)'NNET1=',nnet1
1672 cfmc      call readi(mcmcard,'NNET2',nnet2,4)
1673 cfmc      write(iout,*)'NNET2=',nnet2
1674 cfmc      call readi(mcmcard,'NNET3',nnet3,4)
1675 cfmc      write(iout,*)'NNET3=',nnet3
1676 cfmc      call readi(mcmcard,'ILASTT',ilastt,0)
1677 cfmc      write(iout,*)'ILASTT=',ilastt
1678 cfmc      call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1679 cfmc      write(iout,*)'MAXSTR=',maxstr
1680 cfmc      maxstr_f=maxstr/maxfam
1681 cfmc      write(iout,*)'MAXSTR_F=',maxstr_f
1682 cfmc      call readi(mcmcard,'NMCMF',nmcmf,10)
1683 cfmc      write(iout,*)'NMCMF=',nmcmf
1684 cfmc      call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1685 cfmc      write(iout,*)'IFOCUS=',ifocus
1686 cfmc      call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1687 cfmc      write(iout,*)'NLOCMCMF=',nlocmcmf
1688 cfmc      call readi(mcmcard,'INTPRT',intprt,1000)
1689 cfmc      write(iout,*)'INTPRT=',intprt
1690 cfmc      call readi(mcmcard,'IPRT',iprt,100)
1691 cfmc      write(iout,*)'IPRT=',iprt
1692 cfmc      call readi(mcmcard,'IMAXTR',imaxtr,100)
1693 cfmc      write(iout,*)'IMAXTR=',imaxtr
1694 cfmc      call readi(mcmcard,'MAXEVEN',maxeven,1000)
1695 cfmc      write(iout,*)'MAXEVEN=',maxeven
1696 cfmc      call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1697 cfmc      write(iout,*)'MAXEVEN1=',maxeven1
1698 cfmc      call readi(mcmcard,'INIMIN',inimin,200)
1699 cfmc      write(iout,*)'INIMIN=',inimin
1700 cfmc      call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1701 cfmc      write(iout,*)'NSTEPMCMF=',nstepmcmf
1702 cfmc      call readi(mcmcard,'NTHREAD',nthread,5)
1703 cfmc      write(iout,*)'NTHREAD=',nthread
1704 cfmc      call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1705 cfmc      write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1706 cfmc      call readi(mcmcard,'MAXPERT',maxpert,9)
1707 cfmc      write(iout,*)'MAXPERT=',maxpert
1708 cfmc      call readi(mcmcard,'IRMSD',irmsd,1)
1709 cfmc      write(iout,*)'IRMSD=',irmsd
1710 cfmc      call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1711 cfmc      write(iout,*)'DENEMIN=',denemin
1712 cfmc      call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1713 cfmc      write(iout,*)'RCUT1S=',rcut1s
1714 cfmc      call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1715 cfmc      write(iout,*)'RCUT1E=',rcut1e
1716 cfmc      call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1717 cfmc      write(iout,*)'RCUT2S=',rcut2s
1718 cfmc      call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1719 cfmc      write(iout,*)'RCUT2E=',rcut2e
1720 cfmc      call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1721 cfmc      write(iout,*)'DPERT1=',d_pert1
1722 cfmc      call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1723 cfmc      write(iout,*)'DPERT1A=',d_pert1a
1724 cfmc      call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1725 cfmc      write(iout,*)'DPERT2=',d_pert2
1726 cfmc      call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1727 cfmc      write(iout,*)'DPERT2A=',d_pert2a
1728 cfmc      call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1729 cfmc      write(iout,*)'DPERT2B=',d_pert2b
1730 cfmc      call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1731 cfmc      write(iout,*)'DPERT2C=',d_pert2c
1732 cfmc      d_pert1=deg2rad*d_pert1
1733 cfmc      d_pert1a=deg2rad*d_pert1a
1734 cfmc      d_pert2=deg2rad*d_pert2
1735 cfmc      d_pert2a=deg2rad*d_pert2a
1736 cfmc      d_pert2b=deg2rad*d_pert2b
1737 cfmc      d_pert2c=deg2rad*d_pert2c
1738 cfmc      call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1739 cfmc      write(iout,*)'KT_MCMF1=',kt_mcmf1
1740 cfmc      call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1741 cfmc      write(iout,*)'KT_MCMF2=',kt_mcmf2
1742 cfmc      call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1743 cfmc      write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1744 cfmc      call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1745 cfmc      write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1746 cfmc      call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1747 cfmc      write(iout,*)'RCUTINI=',rcutini
1748 cfmc      call reada(mcmcard,'GRAT',grat,0.5D0)
1749 cfmc      write(iout,*)'GRAT=',grat
1750 cfmc      call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1751 cfmc      write(iout,*)'BIAS_MCMF=',bias_mcmf
1752 cfmc
1753 cfmc      return
1754 cfmc      end
1755 c----------------------------------------------------------------------------
1756       subroutine mcmread
1757       implicit real*8 (a-h,o-z)
1758       include 'DIMENSIONS'
1759       include 'COMMON.MCM'
1760       include 'COMMON.MCE'
1761       include 'COMMON.IOUNITS'
1762       character*80 ucase
1763       character*320 mcmcard
1764       call card_concat(mcmcard)
1765       call readi(mcmcard,'MAXACC',maxacc,100)
1766       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1767       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1768       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1769       call readi(mcmcard,'MAXREPM',maxrepm,200)
1770       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1771       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1772       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1773       call reada(mcmcard,'E_UP',e_up,5.0D0)
1774       call reada(mcmcard,'DELTE',delte,0.1D0)
1775       call readi(mcmcard,'NSWEEP',nsweep,5)
1776       call readi(mcmcard,'NSTEPH',nsteph,0)
1777       call readi(mcmcard,'NSTEPC',nstepc,0)
1778       call reada(mcmcard,'TMIN',tmin,298.0D0)
1779       call reada(mcmcard,'TMAX',tmax,298.0D0)
1780       call readi(mcmcard,'NWINDOW',nwindow,0)
1781       call readi(mcmcard,'PRINT_MC',print_mc,0)
1782       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1783       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1784       ent_read=(index(mcmcard,'ENT_READ').gt.0)
1785       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1786       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1787       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1788       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1789       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1790       if (nwindow.gt.0) then
1791         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1792         do i=1,nwindow
1793           winlen(i)=winend(i)-winstart(i)+1
1794         enddo
1795       endif
1796       if (tmax.lt.tmin) tmax=tmin
1797       if (tmax.eq.tmin) then
1798         nstepc=0
1799         nsteph=0
1800       endif
1801       if (nstepc.gt.0 .and. nsteph.gt.0) then
1802         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1803         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1804       endif
1805 C Probabilities of different move types
1806       sumpro_type(0)=0.0D0
1807       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1808       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1809       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1810       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
1811       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1812       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1813       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1814       do i=1,MaxMoveType
1815         print *,'i',i,' sumprotype',sumpro_type(i)
1816         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1817         print *,'i',i,' sumprotype',sumpro_type(i)
1818       enddo
1819       return
1820       end
1821 c----------------------------------------------------------------------------
1822       subroutine read_minim
1823       implicit real*8 (a-h,o-z)
1824       include 'DIMENSIONS'
1825       include 'COMMON.MINIM'
1826       include 'COMMON.IOUNITS'
1827       character*80 ucase
1828       character*320 minimcard
1829       call card_concat(minimcard)
1830       call readi(minimcard,'MAXMIN',maxmin,2000)
1831       call readi(minimcard,'MAXFUN',maxfun,5000)
1832       call readi(minimcard,'MINMIN',minmin,maxmin)
1833       call readi(minimcard,'MINFUN',minfun,maxmin)
1834       call reada(minimcard,'TOLF',tolf,1.0D-2)
1835       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1836       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1837       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1838       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1839       write (iout,'(/80(1h*)/20x,a/80(1h*))')
1840      &         'Options in energy minimization:'
1841       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1842      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1843      & 'MinMin:',MinMin,' MinFun:',MinFun,
1844      & ' TolF:',TolF,' RTolF:',RTolF
1845       return
1846       end
1847 c----------------------------------------------------------------------------
1848       subroutine read_angles(kanal,*)
1849       implicit real*8 (a-h,o-z)
1850       include 'DIMENSIONS'
1851       include 'COMMON.GEO'
1852       include 'COMMON.VAR'
1853       include 'COMMON.CHAIN'
1854       include 'COMMON.IOUNITS'
1855       include 'COMMON.CONTROL'
1856 c Read angles from input
1857 c
1858        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1859        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1860        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1861        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1862
1863        do i=1,nres
1864 c 9/7/01 avoid 180 deg valence angle
1865         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1866 c
1867         theta(i)=deg2rad*theta(i)
1868         phi(i)=deg2rad*phi(i)
1869         alph(i)=deg2rad*alph(i)
1870         omeg(i)=deg2rad*omeg(i)
1871        enddo
1872       return
1873    10 return1
1874       end
1875 c----------------------------------------------------------------------------
1876       subroutine reada(rekord,lancuch,wartosc,default)
1877       implicit none
1878       character*(*) rekord,lancuch
1879       double precision wartosc,default
1880       integer ilen,iread
1881       external ilen
1882       iread=index(rekord,lancuch)
1883       if (iread.eq.0) then
1884         wartosc=default
1885         return
1886       endif
1887       iread=iread+ilen(lancuch)+1
1888       read (rekord(iread:),*,err=10,end=10) wartosc
1889       return
1890   10  wartosc=default
1891       return
1892       end
1893 c----------------------------------------------------------------------------
1894       subroutine readi(rekord,lancuch,wartosc,default)
1895       implicit none
1896       character*(*) rekord,lancuch
1897       integer wartosc,default
1898       integer ilen,iread
1899       external ilen
1900       iread=index(rekord,lancuch)
1901       if (iread.eq.0) then
1902         wartosc=default
1903         return
1904       endif
1905       iread=iread+ilen(lancuch)+1
1906       read (rekord(iread:),*,err=10,end=10) wartosc
1907       return
1908   10  wartosc=default
1909       return
1910       end
1911 c----------------------------------------------------------------------------
1912       subroutine multreadi(rekord,lancuch,tablica,dim,default)
1913       implicit none
1914       integer dim,i
1915       integer tablica(dim),default
1916       character*(*) rekord,lancuch
1917       character*80 aux
1918       integer ilen,iread
1919       external ilen
1920       do i=1,dim
1921         tablica(i)=default
1922       enddo
1923       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1924       if (iread.eq.0) return
1925       iread=iread+ilen(lancuch)+1
1926       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1927    10 return
1928       end
1929 c----------------------------------------------------------------------------
1930       subroutine multreada(rekord,lancuch,tablica,dim,default)
1931       implicit none
1932       integer dim,i
1933       double precision tablica(dim),default
1934       character*(*) rekord,lancuch
1935       character*80 aux
1936       integer ilen,iread
1937       external ilen
1938       do i=1,dim
1939         tablica(i)=default
1940       enddo
1941       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1942       if (iread.eq.0) return
1943       iread=iread+ilen(lancuch)+1
1944       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1945    10 return
1946       end
1947 c----------------------------------------------------------------------------
1948       subroutine openunits
1949       implicit real*8 (a-h,o-z)
1950       include 'DIMENSIONS'
1951 #ifdef MPI
1952       include 'mpif.h'
1953       character*16 form,nodename
1954       integer nodelen
1955 #endif
1956       include 'COMMON.SETUP'
1957       include 'COMMON.IOUNITS'
1958       include 'COMMON.MD'
1959       include 'COMMON.CONTROL'
1960       integer lenpre,lenpot,ilen,lentmp
1961       external ilen
1962       character*3 out1file_text,ucase
1963       character*3 ll
1964       external ucase
1965       tmpdir=""
1966       out1file_text="YES"
1967 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1968       call getenv_loc("PREFIX",prefix)
1969       pref_orig = prefix
1970       call getenv_loc("POT",pot)
1971       call getenv_loc("DIRTMP",tmpdir)
1972       call getenv_loc("CURDIR",curdir)
1973       call getenv_loc("OUT1FILE",out1file_text)
1974 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1975       out1file_text=ucase(out1file_text)
1976       if (out1file_text(1:1).eq."Y") then
1977         out1file=.true.
1978       else
1979         out1file=fg_rank.gt.0
1980       endif
1981       lenpre=ilen(prefix)
1982       lenpot=ilen(pot)
1983       lentmp=ilen(tmpdir)
1984       if (lentmp.gt.0) then
1985           write (*,'(80(1h!))')
1986           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
1987           write (*,'(80(1h!))')
1988           write (*,*)"All output files will be on node /tmp directory."
1989 #ifdef MPI
1990         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1991         if (me.eq.king) then
1992           write (*,*) "The master node is ",nodename
1993         else if (fg_rank.eq.0) then
1994           write (*,*) "I am the CG slave node ",nodename
1995         else
1996           write (*,*) "I am the FG slave node ",nodename
1997         endif
1998 #endif
1999         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2000         lenpre = lentmp+lenpre+1
2001       endif
2002       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2003 C Get the names and open the input files
2004 #if defined(WINIFL) || defined(WINPGI)
2005       open(1,file=pref_orig(:ilen(pref_orig))//
2006      &  '.inp',status='old',readonly,shared)
2007        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2008 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2009 C Get parameter filenames and open the parameter files.
2010       call getenv_loc('BONDPAR',bondname)
2011       open (ibond,file=bondname,status='old',readonly,shared)
2012       call getenv_loc('THETPAR',thetname)
2013       open (ithep,file=thetname,status='old',readonly,shared)
2014       call getenv_loc('ROTPAR',rotname)
2015       open (irotam,file=rotname,status='old',readonly,shared)
2016       call getenv_loc('TORPAR',torname)
2017       open (itorp,file=torname,status='old',readonly,shared)
2018       call getenv_loc('TORDPAR',tordname)
2019       open (itordp,file=tordname,status='old',readonly,shared)
2020       call getenv_loc('FOURIER',fouriername)
2021       open (ifourier,file=fouriername,status='old',readonly,shared)
2022       call getenv_loc('ELEPAR',elename)
2023       open (ielep,file=elename,status='old',readonly,shared)
2024       call getenv_loc('SIDEPAR',sidename)
2025       open (isidep,file=sidename,status='old',readonly,shared)
2026       call getenv_loc('LIPTRANPAR',liptranname)
2027       open (iliptranpar,file=liptranname,status='old',readonly,shared)
2028 #elif (defined CRAY) || (defined AIX)
2029       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2030      &  action='read')
2031 c      print *,"Processor",myrank," opened file 1"
2032       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2033 c      print *,"Processor",myrank," opened file 9"
2034 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2035 C Get parameter filenames and open the parameter files.
2036       call getenv_loc('BONDPAR',bondname)
2037       open (ibond,file=bondname,status='old',action='read')
2038 c      print *,"Processor",myrank," opened file IBOND"
2039       call getenv_loc('THETPAR',thetname)
2040       open (ithep,file=thetname,status='old',action='read')
2041 #ifndef CRYST_THETA
2042       call getenv_loc('THETPARPDB',thetname_pdb)
2043       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2044 #endif
2045 c      print *,"Processor",myrank," opened file ITHEP"
2046       call getenv_loc('ROTPAR',rotname)
2047       open (irotam,file=rotname,status='old',action='read')
2048 #ifndef CRYST_SC
2049       call getenv_loc('ROTPARPDB',rotname_pdb)
2050       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2051 #endif
2052 c      print *,"Processor",myrank," opened file IROTAM"
2053       call getenv_loc('TORPAR',torname)
2054       open (itorp,file=torname,status='old',action='read')
2055 c      print *,"Processor",myrank," opened file ITORP"
2056       call getenv_loc('TORDPAR',tordname)
2057       open (itordp,file=tordname,status='old',action='read')
2058 c      print *,"Processor",myrank," opened file ITORDP"
2059       call getenv_loc('SCCORPAR',sccorname)
2060       open (isccor,file=sccorname,status='old',action='read')
2061 c      print *,"Processor",myrank," opened file ISCCOR"
2062       call getenv_loc('FOURIER',fouriername)
2063       open (ifourier,file=fouriername,status='old',action='read')
2064 c      print *,"Processor",myrank," opened file IFOURIER"
2065       call getenv_loc('ELEPAR',elename)
2066       open (ielep,file=elename,status='old',action='read')
2067 c      print *,"Processor",myrank," opened file IELEP"
2068       call getenv_loc('SIDEPAR',sidename)
2069       open (isidep,file=sidename,status='old',action='read')
2070       call getenv_loc('LIPTRANPAR',liptranname)
2071       open (iliptranpar,file=liptranname,status='old',action='read')
2072 c      print *,"Processor",myrank," opened file ISIDEP"
2073 c      print *,"Processor",myrank," opened parameter files"
2074 #elif (defined G77)
2075       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2076       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2077 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2078 C Get parameter filenames and open the parameter files.
2079       call getenv_loc('BONDPAR',bondname)
2080       open (ibond,file=bondname,status='old')
2081       call getenv_loc('THETPAR',thetname)
2082       open (ithep,file=thetname,status='old')
2083 #ifndef CRYST_THETA
2084       call getenv_loc('THETPARPDB',thetname_pdb)
2085       open (ithep_pdb,file=thetname_pdb,status='old')
2086 #endif
2087       call getenv_loc('ROTPAR',rotname)
2088       open (irotam,file=rotname,status='old')
2089 #ifndef CRYST_SC
2090       call getenv_loc('ROTPARPDB',rotname_pdb)
2091       open (irotam_pdb,file=rotname_pdb,status='old')
2092 #endif
2093       call getenv_loc('TORPAR',torname)
2094       open (itorp,file=torname,status='old')
2095       call getenv_loc('TORDPAR',tordname)
2096       open (itordp,file=tordname,status='old')
2097       call getenv_loc('SCCORPAR',sccorname)
2098       open (isccor,file=sccorname,status='old')
2099       call getenv_loc('FOURIER',fouriername)
2100       open (ifourier,file=fouriername,status='old')
2101       call getenv_loc('ELEPAR',elename)
2102       open (ielep,file=elename,status='old')
2103       call getenv_loc('SIDEPAR',sidename)
2104       open (isidep,file=sidename,status='old')
2105       call getenv_loc('LIPTRANPAR',liptranname)
2106       open (iliptranpar,file=liptranname,status='old')
2107 #else
2108       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2109      &  readonly)
2110        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2111 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2112 C Get parameter filenames and open the parameter files.
2113       call getenv_loc('BONDPAR',bondname)
2114       open (ibond,file=bondname,status='old',readonly)
2115       call getenv_loc('THETPAR',thetname)
2116       open (ithep,file=thetname,status='old',readonly)
2117       call getenv_loc('ROTPAR',rotname)
2118       open (irotam,file=rotname,status='old',readonly)
2119       call getenv_loc('TORPAR',torname)
2120       open (itorp,file=torname,status='old',readonly)
2121       call getenv_loc('TORDPAR',tordname)
2122       open (itordp,file=tordname,status='old',readonly)
2123       call getenv_loc('SCCORPAR',sccorname)
2124       open (isccor,file=sccorname,status='old',readonly)
2125 #ifndef CRYST_THETA
2126       call getenv_loc('THETPARPDB',thetname_pdb)
2127       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2128 #endif
2129       call getenv_loc('FOURIER',fouriername)
2130       open (ifourier,file=fouriername,status='old',readonly)
2131       call getenv_loc('ELEPAR',elename)
2132       open (ielep,file=elename,status='old',readonly)
2133       call getenv_loc('SIDEPAR',sidename)
2134       open (isidep,file=sidename,status='old',readonly)
2135       call getenv_loc('LIPTRANPAR',liptranname)
2136       open (iliptranpar,file=liptranname,status='old',action='read')
2137 #ifndef CRYST_SC
2138       call getenv_loc('ROTPARPDB',rotname_pdb)
2139       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2140 #endif
2141 #endif
2142 #ifdef TUBE
2143       call getenv_loc('TUBEPAR',tubename)
2144 #if defined(WINIFL) || defined(WINPGI)
2145       open (itube,file=tubename,status='old',readonly,shared)
2146 #elif (defined CRAY)  || (defined AIX)
2147       open (itube,file=tubename,status='old',action='read')
2148 #elif (defined G77)
2149       open (itube,file=tubename,status='old')
2150 #else
2151       open (itube,file=tubename,status='old',readonly)
2152 #endif
2153 #endif
2154 #ifndef OLDSCP
2155 C
2156 C 8/9/01 In the newest version SCp interaction constants are read from a file
2157 C Use -DOLDSCP to use hard-coded constants instead.
2158 C
2159       call getenv_loc('SCPPAR',scpname)
2160 #if defined(WINIFL) || defined(WINPGI)
2161       open (iscpp,file=scpname,status='old',readonly,shared)
2162 #elif (defined CRAY)  || (defined AIX)
2163       open (iscpp,file=scpname,status='old',action='read')
2164 #elif (defined G77)
2165       open (iscpp,file=scpname,status='old')
2166 #else
2167       open (iscpp,file=scpname,status='old',readonly)
2168 #endif
2169 #endif
2170       call getenv_loc('PATTERN',patname)
2171 #if defined(WINIFL) || defined(WINPGI)
2172       open (icbase,file=patname,status='old',readonly,shared)
2173 #elif (defined CRAY)  || (defined AIX)
2174       open (icbase,file=patname,status='old',action='read')
2175 #elif (defined G77)
2176       open (icbase,file=patname,status='old')
2177 #else
2178       open (icbase,file=patname,status='old',readonly)
2179 #endif
2180 #ifdef MPI
2181 C Open output file only for CG processes
2182 c      print *,"Processor",myrank," fg_rank",fg_rank
2183       if (fg_rank.eq.0) then
2184
2185       if (nodes.eq.1) then
2186         npos=3
2187       else
2188         npos = dlog10(dfloat(nodes-1))+1
2189       endif
2190       if (npos.lt.3) npos=3
2191       write (liczba,'(i1)') npos
2192       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2193      &  //')'
2194       write (liczba,form) me
2195       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2196      &  liczba(:ilen(liczba))
2197       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2198      &  //'.int'
2199       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2200      &  //'.pdb'
2201       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2202      &  liczba(:ilen(liczba))//'.mol2'
2203       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2204      &  liczba(:ilen(liczba))//'.stat'
2205       if (lentmp.gt.0)
2206      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2207      &      //liczba(:ilen(liczba))//'.stat')
2208       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2209      &  //'.rst'
2210       if(usampl) then
2211           qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2212      & liczba(:ilen(liczba))//'.const'
2213       endif
2214
2215       endif
2216 #else
2217       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2218       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2219       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2220       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2221       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2222       if (lentmp.gt.0)
2223      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2224      &    //'.stat')
2225       rest2name=prefix(:ilen(prefix))//'.rst'
2226       if(usampl) then
2227          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2228       endif
2229 #endif
2230 #if defined(AIX) || defined(PGI) || defined(CRAY)
2231       if (me.eq.king .or. .not. out1file) then
2232          open(iout,file=outname,status='unknown')
2233       else
2234          open(iout,file="/dev/null",status="unknown")
2235       endif
2236 c#define DEBUG
2237 #ifdef DEBUG
2238       if (fg_rank.gt.0) then
2239         write (liczba,'(i3.3)') myrank/nfgtasks
2240         write (ll,'(bz,i3.3)') fg_rank
2241         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2242      &   status='unknown')
2243       endif
2244 #endif
2245 c#undef DEBUG
2246       if(me.eq.king) then
2247        open(igeom,file=intname,status='unknown',position='append')
2248        open(ipdb,file=pdbname,status='unknown')
2249        open(imol2,file=mol2name,status='unknown')
2250        open(istat,file=statname,status='unknown',position='append')
2251       else
2252 c1out       open(iout,file=outname,status='unknown')
2253       endif
2254 #else
2255       if (me.eq.king .or. .not.out1file)
2256      &    open(iout,file=outname,status='unknown')
2257 #ifdef DEBUG
2258       if (fg_rank.gt.0) then
2259         write (liczba,'(i3.3)') myrank/nfgtasks
2260         write (ll,'(bz,i3.3)') fg_rank
2261         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2262      &   status='unknown')
2263       endif
2264 #endif
2265       if(me.eq.king) then
2266        open(igeom,file=intname,status='unknown',access='append')
2267        open(ipdb,file=pdbname,status='unknown')
2268        open(imol2,file=mol2name,status='unknown')
2269        open(istat,file=statname,status='unknown',access='append')
2270       else
2271 c1out       open(iout,file=outname,status='unknown')
2272       endif
2273 #endif
2274       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2275       csa_seed=prefix(:lenpre)//'.CSA.seed'
2276       csa_history=prefix(:lenpre)//'.CSA.history'
2277       csa_bank=prefix(:lenpre)//'.CSA.bank'
2278       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2279       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2280       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2281 c!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2282       csa_int=prefix(:lenpre)//'.int'
2283       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2284       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2285       csa_in=prefix(:lenpre)//'.CSA.in'
2286 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2287 C Write file names
2288       if (me.eq.king)then
2289       write (iout,'(80(1h-))')
2290       write (iout,'(30x,a)') "FILE ASSIGNMENT"
2291       write (iout,'(80(1h-))')
2292       write (iout,*) "Input file                      : ",
2293      &  pref_orig(:ilen(pref_orig))//'.inp'
2294       write (iout,*) "Output file                     : ",
2295      &  outname(:ilen(outname))
2296       write (iout,*)
2297       write (iout,*) "Sidechain potential file        : ",
2298      &  sidename(:ilen(sidename))
2299 #ifndef OLDSCP
2300       write (iout,*) "SCp potential file              : ",
2301      &  scpname(:ilen(scpname))
2302 #endif
2303       write (iout,*) "Electrostatic potential file    : ",
2304      &  elename(:ilen(elename))
2305       write (iout,*) "Cumulant coefficient file       : ",
2306      &  fouriername(:ilen(fouriername))
2307       write (iout,*) "Torsional parameter file        : ",
2308      &  torname(:ilen(torname))
2309       write (iout,*) "Double torsional parameter file : ",
2310      &  tordname(:ilen(tordname))
2311       write (iout,*) "SCCOR parameter file : ",
2312      &  sccorname(:ilen(sccorname))
2313       write (iout,*) "Bond & inertia constant file    : ",
2314      &  bondname(:ilen(bondname))
2315       write (iout,*) "Bending-torsion parameter file  : ",
2316      &  thetname(:ilen(thetname))
2317       write (iout,*) "Rotamer parameter file          : ",
2318      &  rotname(:ilen(rotname))
2319       write (iout,*) "Thetpdb parameter file          : ",
2320      &  thetname_pdb(:ilen(thetname_pdb))
2321       write (iout,*) "Threading database              : ",
2322      &  patname(:ilen(patname))
2323       if (lentmp.ne.0)
2324      &write (iout,*)" DIRTMP                          : ",
2325      &  tmpdir(:lentmp)
2326       write (iout,'(80(1h-))')
2327       endif
2328       return
2329       end
2330 c----------------------------------------------------------------------------
2331       subroutine card_concat(card)
2332       implicit real*8 (a-h,o-z)
2333       include 'DIMENSIONS'
2334       include 'COMMON.IOUNITS'
2335       character*(*) card
2336       character*80 karta,ucase
2337       external ilen
2338       read (inp,'(a)') karta
2339       karta=ucase(karta)
2340       card=' '
2341       do while (karta(80:80).eq.'&')
2342         card=card(:ilen(card)+1)//karta(:79)
2343         read (inp,'(a)') karta
2344         karta=ucase(karta)
2345       enddo
2346       card=card(:ilen(card)+1)//karta
2347       return
2348       end
2349 c----------------------------------------------------------------------------------
2350       subroutine readrst
2351       implicit real*8 (a-h,o-z)
2352       include 'DIMENSIONS'
2353       include 'COMMON.CHAIN'
2354       include 'COMMON.IOUNITS'
2355       include 'COMMON.MD'
2356       open(irest2,file=rest2name,status='unknown')
2357       read(irest2,*) totT,EK,potE,totE,t_bath
2358       totTafm=totT
2359       do i=1,2*nres
2360          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2361       enddo
2362       do i=1,2*nres
2363          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2364       enddo
2365       if(usampl) then
2366              read (irest2,*) iset
2367       endif
2368       close(irest2)
2369       return
2370       end
2371 c---------------------------------------------------------------------------------
2372       subroutine read_fragments
2373       implicit real*8 (a-h,o-z)
2374       include 'DIMENSIONS'
2375 #ifdef MPI
2376       include 'mpif.h'
2377 #endif
2378       include 'COMMON.SETUP'
2379       include 'COMMON.CHAIN'
2380       include 'COMMON.IOUNITS'
2381       include 'COMMON.MD'
2382       include 'COMMON.CONTROL'
2383       read(inp,*) nset,nfrag,npair,nfrag_back
2384       loc_qlike=(nfrag_back.lt.0)
2385       nfrag_back=iabs(nfrag_back)
2386       if(me.eq.king.or..not.out1file)
2387      & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2388      &  " nfrag_back",nfrag_back," loc_qlike",loc_qlike
2389       do iset=1,nset
2390          read(inp,*) mset(iset)
2391        do i=1,nfrag
2392          read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2393      &     qinfrag(i,iset)
2394          if(me.eq.king.or..not.out1file)
2395      &    write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2396      &     ifrag(2,i,iset), qinfrag(i,iset)
2397        enddo
2398        do i=1,npair
2399         read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2400      &    qinpair(i,iset)
2401         if(me.eq.king.or..not.out1file)
2402      &   write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2403      &    ipair(2,i,iset), qinpair(i,iset)
2404        enddo
2405        if (loc_qlike) then
2406        do i=1,nfrag_back
2407         read(inp,*) wfrag_back(1,i,iset),qin_back(1,i,iset),
2408      &     wfrag_back(2,i,iset),qin_back(2,i,iset),
2409      &     wfrag_back(3,i,iset),qin_back(3,i,iset),
2410      &     ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2411         if(me.eq.king.or..not.out1file)
2412      &   write(iout,*) "A",i,wfrag_back(1,i,iset),qin_back(2,i,iset),
2413      &   wfrag_back(2,i,iset),qin_back(3,i,iset),
2414      &   wfrag_back(3,i,iset),qin_back(3,i,iset),
2415      &   ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2416        enddo
2417        else
2418        do i=1,nfrag_back
2419         read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2420      &     wfrag_back(3,i,iset),
2421      &     ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2422         if(me.eq.king.or..not.out1file)
2423      &   write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2424      &   wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2425        enddo
2426        endif
2427       enddo
2428       return
2429       end
2430 C---------------------------------------------------------------------------
2431       subroutine read_afminp
2432             implicit real*8 (a-h,o-z)
2433       include 'DIMENSIONS'
2434 #ifdef MPI
2435       include 'mpif.h'
2436 #endif
2437       include 'COMMON.SETUP'
2438       include 'COMMON.CONTROL'
2439       include 'COMMON.CHAIN'
2440       include 'COMMON.IOUNITS'
2441       include 'COMMON.SBRIDGE'
2442       character*320 afmcard
2443       print *, "wchodze"
2444       call card_concat(afmcard)
2445       call readi(afmcard,"BEG",afmbeg,0)
2446       call readi(afmcard,"END",afmend,0)
2447       call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2448       call reada(afmcard,"VEL",velAFMconst,0.0d0)
2449       print *,'FORCE=' ,forceAFMconst
2450 CCCC NOW PROPERTIES FOR AFM
2451        distafminit=0.0d0
2452        do i=1,3
2453         distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2454        enddo
2455         distafminit=dsqrt(distafminit)
2456         print *,'initdist',distafminit
2457       return
2458       end
2459 c-------------------------------------------------------------------------------
2460       subroutine read_saxs_constr
2461       implicit real*8 (a-h,o-z)
2462       include 'DIMENSIONS'
2463 #ifdef MPI
2464       include 'mpif.h'
2465 #endif
2466       include 'COMMON.SETUP'
2467       include 'COMMON.CONTROL'
2468       include 'COMMON.CHAIN'
2469       include 'COMMON.IOUNITS'
2470       include 'COMMON.SBRIDGE'
2471       double precision cm(3)
2472 c      read(inp,*) nsaxs
2473       write (iout,*) "Calling read_saxs nsaxs",nsaxs
2474       call flush(iout)
2475       if (saxs_mode.eq.0) then
2476 c SAXS distance distribution
2477       do i=1,nsaxs
2478         read(inp,*) distsaxs(i),Psaxs(i)
2479       enddo
2480       Cnorm = 0.0d0
2481       do i=1,nsaxs
2482         Cnorm = Cnorm + Psaxs(i)
2483       enddo
2484       write (iout,*) "Cnorm",Cnorm
2485       do i=1,nsaxs
2486         Psaxs(i)=Psaxs(i)/Cnorm
2487       enddo
2488       write (iout,*) "Normalized distance distribution from SAXS"
2489       do i=1,nsaxs
2490         write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
2491       enddo
2492       Wsaxs0=0.0d0
2493       do i=1,nsaxs
2494         Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
2495       enddo
2496       write (iout,*) "Wsaxs0",Wsaxs0
2497       else
2498 c SAXS "spheres".
2499       do i=1,nsaxs
2500         read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
2501       enddo
2502       do j=1,3
2503         cm(j)=0.0d0
2504       enddo
2505       do i=1,nsaxs
2506         do j=1,3
2507           cm(j)=cm(j)+Csaxs(j,i)
2508         enddo
2509       enddo
2510       do j=1,3
2511         cm(j)=cm(j)/nsaxs
2512       enddo
2513       do i=1,nsaxs
2514         do j=1,3
2515           Csaxs(j,i)=Csaxs(j,i)-cm(j)
2516         enddo
2517       enddo
2518       write (iout,*) "SAXS sphere coordinates"
2519       do i=1,nsaxs
2520         write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
2521       enddo
2522       endif
2523       return
2524       end
2525
2526 c-------------------------------------------------------------------------------
2527       subroutine read_dist_constr
2528       implicit real*8 (a-h,o-z)
2529       include 'DIMENSIONS'
2530 #ifdef MPI
2531       include 'mpif.h'
2532 #endif
2533       include 'COMMON.SETUP'
2534       include 'COMMON.CONTROL'
2535       include 'COMMON.CHAIN'
2536       include 'COMMON.IOUNITS'
2537       include 'COMMON.SBRIDGE'
2538       include 'COMMON.INTERACT'
2539       integer ifrag_(2,100),ipair_(2,1000)
2540       double precision wfrag_(100),wpair_(1000)
2541       character*5000 controlcard
2542       logical normalize,next
2543       integer restr_type
2544       double precision scal_bfac
2545       double precision xlink(4,0:4) /
2546 c           a          b       c     sigma
2547      &   0.0d0,0.0d0,0.0d0,0.0d0,                             ! default, no xlink potential
2548      &   0.00305218d0,9.46638d0,4.68901d0,4.74347d0,          ! ZL
2549      &   0.00214928d0,12.7517d0,0.00375009d0,6.13477d0,       ! ADH
2550      &   0.00184547d0,11.2678d0,0.00140292d0,7.00868d0,       ! PDH
2551      &   0.000161786d0,6.29273d0,4.40993d0,7.13956d0    /     ! DSS
2552 c      print *, "WCHODZE"
2553       write (iout,*) "Calling read_dist_constr"
2554 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2555 c      call flush(iout)
2556       restr_on_coord=.false.
2557       next=.true.
2558
2559       npeak=0
2560       ipeak=0
2561       nhpb_peak=0
2562
2563       DO WHILE (next)
2564
2565       call card_concat(controlcard)
2566       next = index(controlcard,"NEXT").gt.0
2567       call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
2568       write (iout,*) "restr_type",restr_type
2569       call readi(controlcard,"NFRAG",nfrag_,0)
2570       call readi(controlcard,"NFRAG",nfrag_,0)
2571       call readi(controlcard,"NPAIR",npair_,0)
2572       call readi(controlcard,"NDIST",ndist_,0)
2573       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2574       call reada(controlcard,'SCAL_BFAC',scal_bfac,1.0d0)
2575       if (restr_type.eq.10)
2576      &  call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
2577       if (restr_type.eq.12)
2578      &  call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
2579       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2580       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2581       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2582       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2583       normalize = index(controlcard,"NORMALIZE").gt.0
2584       write (iout,*) "WBOLTZD",wboltzd
2585       write (iout,*) "SCAL_PEAK",scal_peak
2586       write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2587       write (iout,*) "IFRAG"
2588       do i=1,nfrag_
2589         write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2590       enddo
2591       write (iout,*) "IPAIR"
2592       do i=1,npair_
2593         write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2594       enddo
2595       if (nfrag_.gt.0 .or. restr_type.eq.4 .or. restr_type.eq.5)
2596      & write (iout,*)
2597      &   "Distance restraints as generated from reference structure"
2598       do i=1,nfrag_
2599         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2600         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2601      &    ifrag_(2,i)=nstart_sup+nsup-1
2602 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2603 c        call flush(iout)
2604         if (wfrag_(i).eq.0.0d0) cycle
2605         do j=ifrag_(1,i),ifrag_(2,i)-1
2606           do k=j+1,ifrag_(2,i)
2607 c            write (iout,*) "j",j," k",k
2608             ddjk=dist(j,k)
2609             if (restr_type.eq.1) then
2610               nhpb=nhpb+1
2611               irestr_type(nhpb)=1
2612               ihpb(nhpb)=j
2613               jhpb(nhpb)=k
2614               dhpb(nhpb)=ddjk
2615               forcon(nhpb)=wfrag_(i)
2616             else if (constr_dist.eq.2) then
2617               if (ddjk.le.dist_cut) then
2618                 nhpb=nhpb+1
2619                 irestr_type(nhpb)=1
2620                 ihpb(nhpb)=j
2621                 jhpb(nhpb)=k
2622                 dhpb(nhpb)=ddjk
2623                 forcon(nhpb)=wfrag_(i)
2624               endif
2625             else if (restr_type.eq.3) then
2626               nhpb=nhpb+1
2627               irestr_type(nhpb)=1
2628               ihpb(nhpb)=j
2629               jhpb(nhpb)=k
2630               dhpb(nhpb)=ddjk
2631               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2632             endif
2633 #ifdef MPI
2634             if (.not.out1file .or. me.eq.king)
2635      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2636      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2637 #else
2638             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2639      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2640 #endif
2641           enddo
2642         enddo
2643       enddo
2644       do i=1,npair_
2645         if (wpair_(i).eq.0.0d0) cycle
2646         ii = ipair_(1,i)
2647         jj = ipair_(2,i)
2648         if (ii.gt.jj) then
2649           itemp=ii
2650           ii=jj
2651           jj=itemp
2652         endif
2653         do j=ifrag_(1,ii),ifrag_(2,ii)
2654           do k=ifrag_(1,jj),ifrag_(2,jj)
2655             ddjk=dist(j,k)
2656             if (restr_type.eq.1) then
2657               nhpb=nhpb+1
2658               irestr_type(nhpb)=1
2659               ihpb(nhpb)=j
2660               jhpb(nhpb)=k
2661               dhpb(nhpb)=ddjk
2662               forcon(nhpb)=wpair_(i)
2663             else if (constr_dist.eq.2) then
2664               if (ddjk.le.dist_cut) then
2665                 nhpb=nhpb+1
2666                 irestr_type(nhpb)=1
2667                 ihpb(nhpb)=j
2668                 jhpb(nhpb)=k
2669                 dhpb(nhpb)=ddjk
2670                 forcon(nhpb)=wpair_(i)
2671               endif
2672             else if (restr_type.eq.3) then
2673               nhpb=nhpb+1
2674               irestr_type(nhpb)=1
2675               ihpb(nhpb)=j
2676               jhpb(nhpb)=k
2677               dhpb(nhpb)=ddjk
2678               forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2679             endif
2680 #ifdef MPI
2681             if (.not.out1file .or. me.eq.king)
2682      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2683      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2684 #else
2685             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2686      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2687 #endif
2688           enddo
2689         enddo
2690       enddo
2691
2692 c      print *,ndist_
2693       write (iout,*) "Distance restraints as read from input"
2694       do i=1,ndist_
2695         if (restr_type.eq.12) then
2696           read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2697      &    dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2698      &    ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2699      &    fordepth_peak(nhpb_peak+1),npeak
2700 c          write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2701 c     &    dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2702 c     &    ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2703 c     &    fordepth_peak(nhpb_peak+1),npeak
2704           if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
2705      &      fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
2706           nhpb_peak=nhpb_peak+1
2707           irestr_type_peak(nhpb_peak)=12
2708           if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
2709           ipeak(2,npeak)=i
2710 #ifdef MPI
2711           if (.not.out1file .or. me.eq.king)
2712      &    write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2713      &     nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2714      &     ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2715      &     dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2716      &     fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2717 #else
2718           write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2719      &     nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2720      &     ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2721      &     dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2722      &     fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2723 #endif
2724           if (ibecarb_peak(nhpb_peak).eq.3) then
2725             jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2726           else if (ibecarb_peak(nhpb_peak).eq.2) then
2727             ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2728           else if (ibecarb_peak(nhpb_peak).eq.1) then
2729             ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2730             jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2731           endif
2732         else if (restr_type.eq.11) then
2733           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2734      &     dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2735 c        fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2736           if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
2737           nhpb=nhpb+1
2738           irestr_type(nhpb)=11
2739 #ifdef MPI
2740           if (.not.out1file .or. me.eq.king)
2741      &    write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2742      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2743      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2744 #else
2745           write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2746      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2747      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2748 #endif
2749 c          if (ibecarb(nhpb).gt.0) then
2750 c            ihpb(nhpb)=ihpb(nhpb)+nres
2751 c            jhpb(nhpb)=jhpb(nhpb)+nres
2752 c          endif
2753          if (ibecarb(nhpb).eq.3) then
2754             ihpb(nhpb)=ihpb(nhpb)+nres
2755           else if (ibecarb(nhpb).eq.2) then
2756             ihpb(nhpb)=ihpb(nhpb)+nres
2757           else if (ibecarb(nhpb).eq.1) then
2758             ihpb(nhpb)=ihpb(nhpb)+nres
2759             jhpb(nhpb)=jhpb(nhpb)+nres
2760           endif
2761         else if (restr_type.eq.10) then
2762 c Cross-lonk Markov-like potential
2763           call card_concat(controlcard)
2764           call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
2765           call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
2766           ibecarb(nhpb+1)=0
2767           if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
2768           if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
2769           if (index(controlcard,"ZL").gt.0) then
2770             link_type=1
2771           else if (index(controlcard,"ADH").gt.0) then
2772             link_type=2
2773           else if (index(controlcard,"PDH").gt.0) then
2774             link_type=3
2775           else if (index(controlcard,"DSS").gt.0) then
2776             link_type=4
2777           else
2778             link_type=0
2779           endif
2780           call reada(controlcard,"AXLINK",dhpb(nhpb+1),
2781      &       xlink(1,link_type))
2782           call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
2783      &       xlink(2,link_type))
2784           call reada(controlcard,"CXLINK",fordepth(nhpb+1),
2785      &       xlink(3,link_type))
2786           call reada(controlcard,"SIGMA",forcon(nhpb+1),
2787      &       xlink(4,link_type))
2788           call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
2789 c          read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
2790 c     &      dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2791           if (forcon(nhpb+1).le.0.0d0 .or.
2792      &       (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
2793           nhpb=nhpb+1
2794           irestr_type(nhpb)=10
2795           if (ibecarb(nhpb).eq.3) then
2796             jhpb(nhpb)=jhpb(nhpb)+nres
2797           else if (ibecarb(nhpb).eq.2) then
2798             ihpb(nhpb)=ihpb(nhpb)+nres
2799           else if (ibecarb(nhpb).eq.1) then
2800             ihpb(nhpb)=ihpb(nhpb)+nres
2801             jhpb(nhpb)=jhpb(nhpb)+nres
2802           endif
2803 #ifdef MPI
2804           if (.not.out1file .or. me.eq.king)
2805      &    write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2806      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2807      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2808      &     irestr_type(nhpb)
2809 #else
2810           write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2811      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2812      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2813      &     irestr_type(nhpb)
2814 #endif
2815         else
2816 C        print *,"in else"
2817           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2818      &     dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
2819           if (forcon(nhpb+1).gt.0.0d0) then
2820           nhpb=nhpb+1
2821           if (dhpb1(nhpb).eq.0.0d0) then
2822             irestr_type(nhpb)=1
2823           else
2824             irestr_type(nhpb)=2
2825           endif
2826           if (ibecarb(nhpb).eq.3) then
2827             jhpb(nhpb)=jhpb(nhpb)+nres
2828           else if (ibecarb(nhpb).eq.2) then
2829             ihpb(nhpb)=ihpb(nhpb)+nres
2830           else if (ibecarb(nhpb).eq.1) then
2831             ihpb(nhpb)=ihpb(nhpb)+nres
2832             jhpb(nhpb)=jhpb(nhpb)+nres
2833           endif
2834           if (dhpb(nhpb).eq.0.0d0)
2835      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2836         endif
2837 #ifdef MPI
2838           if (.not.out1file .or. me.eq.king)
2839      &    write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2840      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2841 #else
2842           write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2843      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2844 #endif
2845         endif
2846 C        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2847 C        if (forcon(nhpb+1).gt.0.0d0) then
2848 C          nhpb=nhpb+1
2849 C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2850       enddo
2851
2852       if (restr_type.eq.4) then
2853         write (iout,*) "The BFAC array"
2854         do i=nnt,nct
2855           write (iout,'(i5,f10.5)') i,bfac(i)
2856         enddo
2857         do i=nnt,nct
2858           if (itype(i).eq.ntyp1) cycle
2859           do j=nnt,i-1
2860             if (itype(j).eq.ntyp1) cycle
2861             if (itype(i).eq.10) then
2862               iiend=0
2863             else
2864               iiend=1
2865             endif
2866             if (itype(j).eq.10) then
2867               jjend=0
2868             else
2869               jjend=1
2870             endif
2871             kk=0
2872             do ii=0,iiend
2873             do jj=0,jjend
2874             nhpb=nhpb+1
2875             irestr_type(nhpb)=1
2876             forcon(nhpb)=scal_bfac**2/(bfac(i)**2+bfac(j)**2)
2877             irestr_type(nhpb)=1
2878             ibecarb(nhpb)=kk
2879             if (ibecarb(nhpb).gt.0) ibecarb(nhpb)=4-ibecarb(nhpb)
2880             ihpb(nhpb)=i+nres*ii
2881             jhpb(nhpb)=j+nres*jj
2882             dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
2883 #ifdef MPI
2884             if (.not.out1file .or. me.eq.king) then
2885             write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2886      &       nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2887      &       dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2888      &       irestr_type(nhpb)
2889             endif
2890 #else
2891             write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2892      &       nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2893      &       dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2894      &       irestr_type(nhpb)
2895 #endif
2896             kk=kk+1
2897           enddo
2898           enddo
2899           enddo
2900         enddo
2901       endif
2902
2903       if (restr_type.eq.5) then
2904         restr_on_coord=.true.
2905         do i=nnt,nct
2906           if (itype(i).eq.ntyp1) cycle
2907           bfac(i)=(scal_bfac/bfac(i))**2
2908         enddo
2909       endif
2910
2911       ENDDO ! next
2912
2913       fordepthmax=0.0d0
2914       if (normalize) then
2915         do i=nss+1,nhpb
2916           if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
2917      &      fordepthmax=fordepth(i)
2918         enddo
2919         do i=nss+1,nhpb
2920           if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
2921         enddo
2922       endif
2923       return
2924       end
2925 c-------------------------------------------------------------------------------
2926 #ifdef WINIFL
2927       subroutine flush(iu)
2928       return
2929       end
2930 #endif
2931 #ifdef AIX
2932       subroutine flush(iu)
2933       call flush_(iu)
2934       return
2935       end
2936 #endif
2937 c------------------------------------------------------------------------------
2938       subroutine copy_to_tmp(source)
2939       include "DIMENSIONS"
2940       include "COMMON.IOUNITS"
2941       character*(*) source
2942       character* 256 tmpfile
2943       integer ilen
2944       external ilen
2945       logical ex
2946       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2947       inquire(file=tmpfile,exist=ex)
2948       if (ex) then
2949         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2950      &   " to temporary directory..."
2951         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2952         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2953       endif
2954       return
2955       end
2956 c------------------------------------------------------------------------------
2957       subroutine move_from_tmp(source)
2958       include "DIMENSIONS"
2959       include "COMMON.IOUNITS"
2960       character*(*) source
2961       integer ilen
2962       external ilen
2963       write (*,*) "Moving ",source(:ilen(source)),
2964      & " from temporary directory to working directory"
2965       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2966       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2967       return
2968       end
2969 c------------------------------------------------------------------------------
2970       subroutine random_init(seed)
2971 C
2972 C Initialize random number generator
2973 C
2974       implicit real*8 (a-h,o-z)
2975       include 'DIMENSIONS'
2976 #ifdef MPI
2977       include 'mpif.h'
2978       logical OKRandom, prng_restart
2979       real*8  r1
2980       integer iseed_array(4)
2981 #endif
2982       include 'COMMON.IOUNITS'
2983       include 'COMMON.TIME1'
2984       include 'COMMON.THREAD'
2985       include 'COMMON.SBRIDGE'
2986       include 'COMMON.CONTROL'
2987       include 'COMMON.MCM'
2988       include 'COMMON.MAP'
2989       include 'COMMON.HEADER'
2990       include 'COMMON.CSA'
2991       include 'COMMON.CHAIN'
2992       include 'COMMON.MUCA'
2993       include 'COMMON.MD'
2994       include 'COMMON.FFIELD'
2995       include 'COMMON.SETUP'
2996       iseed=-dint(dabs(seed))
2997       if (iseed.eq.0) then
2998         write (iout,'(/80(1h*)/20x,a/80(1h*))')
2999      &    'Random seed undefined. The program will stop.'
3000         write (*,'(/80(1h*)/20x,a/80(1h*))')
3001      &    'Random seed undefined. The program will stop.'
3002 #ifdef MPI
3003         call mpi_finalize(mpi_comm_world,ierr)
3004 #endif
3005         stop 'Bad random seed.'
3006       endif
3007 #ifdef MPI
3008       if (fg_rank.eq.0) then
3009       seed=seed*(me+1)+1
3010 #ifdef AMD64
3011       if(me.eq.king)
3012      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
3013       OKRandom = prng_restart(me,iseed)
3014 #else
3015       do i=1,4
3016        tmp=65536.0d0**(4-i)
3017        iseed_array(i) = dint(seed/tmp)
3018        seed=seed-iseed_array(i)*tmp
3019       enddo
3020       if(me.eq.king)
3021      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3022      &                 (iseed_array(i),i=1,4)
3023       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3024      &                 (iseed_array(i),i=1,4)
3025       OKRandom = prng_restart(me,iseed_array)
3026 #endif
3027       if (OKRandom) then
3028 c        r1 = prng_next(me)
3029         r1=ran_number(0.0D0,1.0D0)
3030         if(me.eq.king)
3031      &   write (iout,*) 'ran_num',r1
3032         if (r1.lt.0.0d0) OKRandom=.false.
3033       endif
3034       if (.not.OKRandom) then
3035         write (iout,*) 'PRNG IS NOT WORKING!!!'
3036         print *,'PRNG IS NOT WORKING!!!'
3037         if (me.eq.0) then
3038          call flush(iout)
3039          call mpi_abort(mpi_comm_world,error_msg,ierr)
3040          stop
3041         else
3042          write (iout,*) 'too many processors for parallel prng'
3043          write (*,*) 'too many processors for parallel prng'
3044          call flush(iout)
3045          stop
3046         endif
3047       endif
3048       endif
3049 #else
3050       call vrndst(iseed)
3051       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
3052 #endif
3053       return
3054       end