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