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