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