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