source/unres/src_CSA/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       logical file_exist
  12 C Read force-field parameters except weights
  13       call parmread
  14 C Read job setup parameters
  15       call read_control
  16 C Read control parameters for energy minimzation if required
  17       if (minim) call read_minim
  18 C Read MCM control parameters if required
  19 c      if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
  20 C Read MD control parameters if reqjuired
  21 c      if (modecalc.eq.12) call read_MDpar
  22 C Read MREMD control parameters if required
  23 c      if (modecalc.eq.14) then
  24 c         call read_MDpar
  25 c         call read_REMDpar
  26 c      endif
  27 C Read MUCA control parameters if required
  28 c      if (lmuca) call read_muca
  29 C Read CSA control parameters if required (from fort.40 if exists
  30 C otherwise from general input file)
  31       if (modecalc.eq.8) then
  32        inquire (file="fort.40",exist=file_exist)
  33        if (.not.file_exist) call csaread
  34       endif
  35 cfmc      if (modecalc.eq.10) call mcmfread
  36 C Read molecule information, molecule geometry, energy-term weights, and
  37 C restraints if requested
  38       call molread
  39
  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 c      include 'COMMON.THREAD'
  71       include 'COMMON.SBRIDGE'
  72       include 'COMMON.CONTROL'
  73       include 'COMMON.MCM'
  74 c      include 'COMMON.MAP'
  75       include 'COMMON.HEADER'
  76 c      include 'COMMON.CSA'
  77       include 'COMMON.CHAIN'
  78 c      include 'COMMON.MUCA'
  79 c      include 'COMMON.MD'
  80       include 'COMMON.FFIELD'
  81       include 'COMMON.SETUP'
  82       COMMON /MACHSW/ KDIAG,ICORFL,IXDR
  83       character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
  84       character*80 ucase
  85       character*320 controlcard
  86
  87       nglob_csa=0
  88       eglob_csa=1d99
  89       nmin_csa=0
  90       read (INP,'(a)') titel
  91       call card_concat(controlcard)
  92 c      out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
  93 c      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
  94       call reada(controlcard,'SEED',seed,0.0D0)
  95       call random_init(seed)
  96 C Set up the time limit (caution! The time must be input in minutes!)
  97       read_cart=index(controlcard,'READ_CART').gt.0
  98       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
  99       call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
 100       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
 101       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
 102       call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
 103       call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
 104       call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
 105       call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
 106       call reada(controlcard,'DRMS',drms,0.1D0)
 107       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 108        write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
 109        write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
 110        write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
 111        write (iout,'(a,f10.1)')'DRMS    = ',drms
 112        write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
 113        write (iout,'(a,f10.1)') 'Time limit (min):',timlim
 114       endif
 115       call readi(controlcard,'NZ_START',nz_start,0)
 116       call readi(controlcard,'NZ_END',nz_end,0)
 117       call readi(controlcard,'IZ_SC',iz_sc,0)
 118       timlim=60.0D0*timlim
 119       safety = 60.0d0*safety
 120       timem=timlim
 121       modecalc=0
 122       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 123       minim=(index(controlcard,'MINIMIZE').gt.0)
 124       dccart=(index(controlcard,'CART').gt.0)
 125       overlapsc=(index(controlcard,'OVERLAP').gt.0)
 126       overlapsc=.not.overlapsc
 127       searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
 128       searchsc=.not.searchsc
 129       sideadd=(index(controlcard,'SIDEADD').gt.0)
 130       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
 131       outpdb=(index(controlcard,'PDBOUT').gt.0)
 132       outmol2=(index(controlcard,'MOL2OUT').gt.0)
 133       pdbref=(index(controlcard,'PDBREF').gt.0)
 134       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
 135       indpdb=index(controlcard,'PDBSTART')
 136       extconf=(index(controlcard,'EXTCONF').gt.0)
 137       call readi(controlcard,'IPRINT',iprint,0)
 138       call readi(controlcard,'MAXGEN',maxgen,10000)
 139       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
 140       call readi(controlcard,"KDIAG",kdiag,0)
 141       call readi(controlcard,"RESCALE_MODE",rescale_mode,0)
 142       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
 143      & write (iout,*) "RESCALE_MODE",rescale_mode
 144       split_ene=index(controlcard,'SPLIT_ENE').gt.0
 145       if (index(controlcard,'REGULAR').gt.0.0D0) then
 146         call reada(controlcard,'WEIDIS',weidis,0.1D0)
 147         modecalc=1
 148         refstr=.true.
 149       endif
 150       if (index(controlcard,'CHECKGRAD').gt.0) then
 151         modecalc=5
 152         if (index(controlcard,'CART').gt.0) then
 153           icheckgrad=1
 154         elseif (index(controlcard,'CARINT').gt.0) then
 155           icheckgrad=2
 156         else
 157           icheckgrad=3
 158         endif
 159       elseif (index(controlcard,'THREAD').gt.0) then
 160         modecalc=2
 161         call readi(controlcard,'THREAD',nthread,0)
 162         if (nthread.gt.0) then
 163           call reada(controlcard,'WEIDIS',weidis,0.1D0)
 164         else
 165           if (fg_rank.eq.0)
 166      &    write (iout,'(a)')'A number has to follow the THREAD keyword.'
 167           stop 'Error termination in Read_Control.'
 168         endif
 169       else if (index(controlcard,'MCMA').gt.0) then
 170         modecalc=3
 171       else if (index(controlcard,'MCEE').gt.0) then
 172         modecalc=6
 173       else if (index(controlcard,'MULTCONF').gt.0) then
 174         modecalc=4
 175       else if (index(controlcard,'MAP').gt.0) then
 176         modecalc=7
 177         call readi(controlcard,'MAP',nmap,0)
 178       else if (index(controlcard,'CSA').gt.0) then
 179         modecalc=8
 180 crc      else if (index(controlcard,'ZSCORE').gt.0) then
 181 crc
 182 crc  ZSCORE is rm from UNRES, modecalc=9 is available
 183 crc
 184 crc        modecalc=9
 185 cfcm      else if (index(controlcard,'MCMF').gt.0) then
 186 cfmc        modecalc=10
 187       else if (index(controlcard,'SOFTREG').gt.0) then
 188         modecalc=11
 189       else if (index(controlcard,'CHECK_BOND').gt.0) then
 190         modecalc=-1
 191       else if (index(controlcard,'TEST').gt.0) then
 192         modecalc=-2
 193       else if (index(controlcard,'MD').gt.0) then
 194         modecalc=12
 195       else if (index(controlcard,'RE ').gt.0) then
 196         modecalc=14
 197       endif
 198
 199       lmuca=index(controlcard,'MUCA').gt.0
 200       call readi(controlcard,'MUCADYN',mucadyn,0)
 201       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
 202       if (lmuca .and. (me.eq.king .or. .not.out1file ))
 203      & then
 204        write (iout,*) 'MUCADYN=',mucadyn
 205        write (iout,*) 'MUCASMOOTH=',muca_smooth
 206       endif
 207
 208       iscode=index(controlcard,'ONE_LETTER')
 209       indphi=index(controlcard,'PHI')
 210       indback=index(controlcard,'BACK')
 211       iranconf=index(controlcard,'RAND_CONF')
 212       i2ndstr=index(controlcard,'USE_SEC_PRED')
 213       gradout=index(controlcard,'GRADOUT').gt.0
 214       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
 215
 216       if(me.eq.king.or..not.out1file)
 217      & write (iout,'(2a)') diagmeth(kdiag),
 218      &  ' routine used to diagonalize matrices.'
 219       return
 220       end
 221 c------------------------------------------------------------------------------
 222       subroutine molread
 223 C
 224 C Read molecular data.
 225 C
 226       implicit real*8 (a-h,o-z)
 227       include 'DIMENSIONS'
 228 #ifdef MPI
 229       include 'mpif.h'
 230       integer error_msg
 231 #endif
 232       include 'COMMON.IOUNITS'
 233       include 'COMMON.GEO'
 234       include 'COMMON.VAR'
 235       include 'COMMON.INTERACT'
 236       include 'COMMON.LOCAL'
 237       include 'COMMON.NAMES'
 238       include 'COMMON.CHAIN'
 239       include 'COMMON.FFIELD'
 240       include 'COMMON.SBRIDGE'
 241       include 'COMMON.HEADER'
 242       include 'COMMON.CONTROL'
 243 c      include 'COMMON.DBASE'
 244 c      include 'COMMON.THREAD'
 245       include 'COMMON.CONTACTS'
 246       include 'COMMON.TORCNSTR'
 247       include 'COMMON.TIME1'
 248       include 'COMMON.BOUNDS'
 249 c      include 'COMMON.MD'
 250 c      include 'COMMON.REMD'
 251       include 'COMMON.SETUP'
 252       character*4 sequence(maxres)
 253       integer rescode
 254       double precision x(maxvar)
 255       character*256 pdbfile
 256       character*320 weightcard
 257       character*80 weightcard_t,ucase
 258       dimension itype_pdb(maxres)
 259       common /pizda/ itype_pdb
 260       logical seq_comp,fail
 261       double precision energia(0:n_ene)
 262       integer ilen
 263       external ilen
 264 C
 265 C Body
 266 C
 267 C Read weights of the subsequent energy terms.
 268
 269        call card_concat(weightcard)
 270        call reada(weightcard,'WLONG',wlong,1.0D0)
 271        call reada(weightcard,'WSC',wsc,wlong)
 272        call reada(weightcard,'WSCP',wscp,wlong)
 273        call reada(weightcard,'WELEC',welec,1.0D0)
 274        call reada(weightcard,'WVDWPP',wvdwpp,welec)
 275        call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 276        call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 277        call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 278        call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 279        call reada(weightcard,'WTURN3',wturn3,1.0D0)
 280        call reada(weightcard,'WTURN4',wturn4,1.0D0)
 281        call reada(weightcard,'WTURN6',wturn6,1.0D0)
 282        call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 283        call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 284        call reada(weightcard,'WBOND',wbond,1.0D0)
 285        call reada(weightcard,'WTOR',wtor,1.0D0)
 286        call reada(weightcard,'WTORD',wtor_d,1.0D0)
 287        call reada(weightcard,'WANG',wang,1.0D0)
 288        call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 289 C     Juyong
 290        call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
 291        call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
 292        call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
 293        call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
 294 C
 295        call reada(weightcard,'SCAL14',scal14,0.4D0)
 296        call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 297        call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 298        call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 299        call reada(weightcard,'TEMP0',temp0,300.0d0)
 300        if (index(weightcard,'SOFT').gt.0) ipot=6
 301 C 12/1/95 Added weight for the multi-body term WCORR
 302        call reada(weightcard,'WCORRH',wcorr,1.0D0)
 303        if (wcorr4.gt.0.0d0) wcorr=wcorr4
 304        weights(1)=wsc
 305        weights(2)=wscp
 306        weights(3)=welec
 307        weights(4)=wcorr
 308        weights(5)=wcorr5
 309        weights(6)=wcorr6
 310        weights(7)=wel_loc
 311        weights(8)=wturn3
 312        weights(9)=wturn4
 313        weights(10)=wturn6
 314        weights(11)=wang
 315        weights(12)=wscloc
 316        weights(13)=wtor
 317        weights(14)=wtor_d
 318        weights(15)=wstrain
 319        weights(16)=wvdwpp
 320        weights(17)=wbond
 321        weights(18)=scal14
 322        weights(21)=wsccor
 323 C     JUYONG
 324        weights(24)=wdfa_dist
 325        weights(25)=wdfa_tor
 326        weights(26)=wdfa_nei
 327        weights(27)=wdfa_beta
 328 C
 329
 330       if(me.eq.king.or..not.out1file)
 331      & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 332      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
 333      &  wturn4,wturn6,
 334      &  wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
 335
 336    10 format (/'Energy-term weights (unscaled):'//
 337      & 'WSCC=   ',f10.6,' (SC-SC)'/
 338      & 'WSCP=   ',f10.6,' (SC-p)'/
 339      & 'WELEC=  ',f10.6,' (p-p electr)'/
 340      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 341      & 'WBOND=  ',f10.6,' (stretching)'/
 342      & 'WANG=   ',f10.6,' (bending)'/
 343      & 'WSCLOC= ',f10.6,' (SC local)'/
 344      & 'WTOR=   ',f10.6,' (torsional)'/
 345      & 'WTORD=  ',f10.6,' (double torsional)'/
 346      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 347      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 348      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 349      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 350      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 351      & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
 352      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 353      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 354      & 'WTURN6= ',f10.6,' (turns, 6th order)'/
 355      & 'WDFA_D= ',f10.6,' (DFA, distance)'   /
 356      & 'WDFA_T= ',f10.6,' (DFA, torsional)'   /
 357      & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)'   /
 358      & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
 359
 360       if(me.eq.king.or..not.out1file)then
 361        if (wcorr4.gt.0.0d0) then
 362         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
 363      &   'between contact pairs of peptide groups'
 364         write (iout,'(2(a,f5.3/))')
 365      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
 366      &  'Range of quenching the correlation terms:',2*delt_corr
 367        else if (wcorr.gt.0.0d0) then
 368         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
 369      &   'between contact pairs of peptide groups'
 370        endif
 371        write (iout,'(a,f8.3)')
 372      &  'Scaling factor of 1,4 SC-p interactions:',scal14
 373        write (iout,'(a,f8.3)')
 374      &  'General scaling factor of SC-p interactions:',scalscp
 375       endif
 376       r0_corr=cutoff_corr-delt_corr
 377       do i=1,20
 378         aad(i,1)=scalscp*aad(i,1)
 379         aad(i,2)=scalscp*aad(i,2)
 380         bad(i,1)=scalscp*bad(i,1)
 381         bad(i,2)=scalscp*bad(i,2)
 382       enddo
 383 c      call rescale_weights(t_bath)
 384       if(me.eq.king.or..not.out1file)
 385      & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 386      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
 387      &  wturn4,wturn6,
 388      &  wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
 389
 390    22 format (/'Energy-term weights (scaled):'//
 391      & 'WSCC=   ',f10.6,' (SC-SC)'/
 392      & 'WSCP=   ',f10.6,' (SC-p)'/
 393      & 'WELEC=  ',f10.6,' (p-p electr)'/
 394      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 395      & 'WBOND=  ',f10.6,' (stretching)'/
 396      & 'WANG=   ',f10.6,' (bending)'/
 397      & 'WSCLOC= ',f10.6,' (SC local)'/
 398      & 'WTOR=   ',f10.6,' (torsional)'/
 399      & 'WTORD=  ',f10.6,' (double torsional)'/
 400      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 401      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 402      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 403      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 404      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 405      & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
 406      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 407      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 408      & 'WTURN6= ',f10.6,' (turns, 6th order)'/
 409      & 'WDFA_D= ',f10.6,' (DFA, distance)'   /
 410      & 'WDFA_T= ',f10.6,' (DFA, torsional)'   /
 411      & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)'   /
 412      & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
 413
 414       if(me.eq.king.or..not.out1file)
 415      & write (iout,*) "Reference temperature for weights calculation:",
 416      &  temp0
 417       call reada(weightcard,"D0CM",d0cm,3.78d0)
 418       call reada(weightcard,"AKCM",akcm,15.1d0)
 419       call reada(weightcard,"AKTH",akth,11.0d0)
 420       call reada(weightcard,"AKCT",akct,12.0d0)
 421       call reada(weightcard,"V1SS",v1ss,-1.08d0)
 422       call reada(weightcard,"V2SS",v2ss,7.61d0)
 423       call reada(weightcard,"V3SS",v3ss,13.7d0)
 424       call reada(weightcard,"EBR",ebr,-5.50D0)
 425       if(me.eq.king.or..not.out1file) then
 426        write (iout,*) "Parameters of the SS-bond potential:"
 427        write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
 428      & " AKCT",akct
 429        write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
 430        write (iout,*) "EBR",ebr
 431        print *,'indpdb=',indpdb,' pdbref=',pdbref
 432       endif
 433       if (indpdb.gt.0 .or. pdbref) then
 434         read(inp,'(a)') pdbfile
 435         if(me.eq.king.or..not.out1file)
 436      &   write (iout,'(2a)') 'PDB data will be read from file ',
 437      &   pdbfile(:ilen(pdbfile))
 438         open(ipdbin,file=pdbfile,status='old',err=33)
 439         goto 34
 440   33    write (iout,'(a)') 'Error opening PDB file.'
 441         stop
 442   34    continue
 443 c        print *,'Begin reading pdb data'
 444         call readpdb
 445 c        print *,'Finished reading pdb data'
 446         if(me.eq.king.or..not.out1file)
 447      &   write (iout,'(a,i3,a,i3)')'nsup=',nsup,
 448      &   ' nstart_sup=',nstart_sup
 449         do i=1,nres
 450           itype_pdb(i)=itype(i)
 451         enddo
 452         close (ipdbin)
 453         nnt=nstart_sup
 454         nct=nstart_sup+nsup-1
 455         call contact(.false.,ncont_ref,icont_ref,co)
 456
 457         if (sideadd) then
 458 C Following 2 lines for diagnostics; comment out if not needed
 459          write (iout,*) "Before sideadd"
 460          call intout
 461          if(me.eq.king.or..not.out1file)
 462      &    write(iout,*)'Adding sidechains'
 463          maxsi=1000
 464          do i=2,nres-1
 465           iti=itype(i)
 466           if (iti.ne.10) then
 467             nsi=0
 468             fail=.true.
 469             do while (fail.and.nsi.le.maxsi)
 470 c              call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
 471               nsi=nsi+1
 472             enddo
 473             if(fail) write(iout,*)'Adding sidechain failed for res ',
 474      &              i,' after ',nsi,' trials'
 475           endif
 476          enddo
 477         endif
 478 C 10/03/12 Adam: Recalculate coordinates with new side chain positions
 479         call chainbuild
 480 C Following 2 lines for diagnostics; comment out if not needed
 481         write (iout,*) "After sideadd"
 482         call intout
 483       endif
 484
 485       if (indpdb.eq.0) then
 486 C Read sequence if not taken from the pdb file.
 487         read (inp,*) nres
 488 c        print *,'nres=',nres
 489         if (iscode.gt.0) then
 490           read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 491         else
 492           read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 493         endif
 494 C Convert sequence to numeric code
 495         do i=1,nres
 496           itype(i)=rescode(i,sequence(i),iscode)
 497         enddo
 498 C Assign initial virtual bond lengths
 499         do i=2,nres
 500           vbld(i)=vbl
 501           vbld_inv(i)=vblinv
 502         enddo
 503         do i=2,nres-1
 504           vbld(i+nres)=dsc(itype(i))
 505           vbld_inv(i+nres)=dsc_inv(itype(i))
 506 c          write (iout,*) "i",i," itype",itype(i),
 507 c     &      " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
 508         enddo
 509       endif
 510 c      print *,nres
 511 c      print '(20i4)',(itype(i),i=1,nres)
 512       do i=1,nres
 513 #ifdef PROCOR
 514         if (itype(i).eq.21 .or. itype(i+1).eq.21) then
 515 #else
 516         if (itype(i).eq.21) then
 517 #endif
 518           itel(i)=0
 519 #ifdef PROCOR
 520         else if (itype(i+1).ne.20) then
 521 #else
 522         else if (itype(i).ne.20) then
 523 #endif
 524           itel(i)=1
 525         else
 526           itel(i)=2
 527         endif
 528       enddo
 529       if(me.eq.king.or..not.out1file)then
 530        write (iout,*) "ITEL"
 531        do i=1,nres-1
 532          write (iout,*) i,itype(i),itel(i)
 533        enddo
 534        print *,'Call Read_Bridge.'
 535       endif
 536       call read_bridge
 537 C 8/13/98 Set limits to generating the dihedral angles
 538       do i=1,nres
 539         phibound(1,i)=-pi
 540         phibound(2,i)=pi
 541       enddo
 542       read (inp,*) ndih_constr
 543       if (ndih_constr.gt.0) then
 544         read (inp,*) ftors
 545         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
 546         if(me.eq.king.or..not.out1file)then
 547          write (iout,*)
 548      &   'There are',ndih_constr,' constraints on phi angles.'
 549          do i=1,ndih_constr
 550           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
 551          enddo
 552         endif
 553         do i=1,ndih_constr
 554           phi0(i)=deg2rad*phi0(i)
 555           drange(i)=deg2rad*drange(i)
 556         enddo
 557         if(me.eq.king.or..not.out1file)
 558      &   write (iout,*) 'FTORS',ftors
 559         do i=1,ndih_constr
 560           ii = idih_constr(i)
 561           phibound(1,ii) = phi0(i)-drange(i)
 562           phibound(2,ii) = phi0(i)+drange(i)
 563         enddo
 564       endif
 565       nnt=1
 566 #ifdef MPI
 567       if (me.eq.king) then
 568 #endif
 569        write (iout,'(a)') 'Boundaries in phi angle sampling:'
 570        do i=1,nres
 571          write (iout,'(a3,i5,2f10.1)')
 572      &   restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
 573        enddo
 574 #ifdef MP
 575       endif
 576 #endif
 577       nct=nres
 578 cd      print *,'NNT=',NNT,' NCT=',NCT
 579       if (itype(1).eq.21) nnt=2
 580       if (itype(nres).eq.21) nct=nct-1
 581
 582 C     Juyong:READ init_vars
 583 C     Initialize variables!
 584 C     Juyong:READ read_info
 585 C     READ fragment information!!
 586 C     both routines should be in dfa.F file!!
 587
 588       if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
 589      &            wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
 590        call init_dfa_vars
 591        print*, 'init_dfa_vars finished!'
 592        call read_dfa_info
 593        print*, 'read_dfa_info finished!'
 594       endif
 595 C
 596 C
 597
 598
 599       if (pdbref) then
 600         if(me.eq.king.or..not.out1file)
 601      &   write (iout,'(a,i3)') 'nsup=',nsup
 602         nstart_seq=nnt
 603         if (nsup.le.(nct-nnt+1)) then
 604           do i=0,nct-nnt+1-nsup
 605             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 606               nstart_seq=nnt+i
 607               goto 111
 608             endif
 609           enddo
 610           write (iout,'(a)')
 611      &            'Error - sequences to be superposed do not match.'
 612           stop
 613         else
 614           do i=0,nsup-(nct-nnt+1)
 615             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 616      &      then
 617               nstart_sup=nstart_sup+i
 618               nsup=nct-nnt+1
 619               goto 111
 620             endif
 621           enddo
 622           write (iout,'(a)')
 623      &            'Error - sequences to be superposed do not match.'
 624         endif
 625   111   continue
 626         if (nsup.eq.0) nsup=nct-nnt
 627         if (nstart_sup.eq.0) nstart_sup=nnt
 628         if (nstart_seq.eq.0) nstart_seq=nnt
 629         if(me.eq.king.or..not.out1file)
 630      &   write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
 631      &                 ' nstart_seq=',nstart_seq
 632       endif
 633 c--- Zscore rms -------
 634       if (nz_start.eq.0) nz_start=nnt
 635       if (nz_end.eq.0 .and. nsup.gt.0) then
 636         nz_end=nnt+nsup-1
 637       else if (nz_end.eq.0) then
 638         nz_end=nct
 639       endif
 640       if(me.eq.king.or..not.out1file)then
 641        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
 642        write (iout,*) 'IZ_SC=',iz_sc
 643       endif
 644 c----------------------
 645       call init_int_table
 646       if (refstr) then
 647         if (.not.pdbref) then
 648           call read_angles(inp,*38)
 649           goto 39
 650    38     write (iout,'(a)') 'Error reading reference structure.'
 651 #ifdef MPI
 652           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
 653           stop 'Error reading reference structure'
 654 #endif
 655    39     call chainbuild
 656           call setup_var
 657 czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
 658           nstart_sup=nnt
 659           nstart_seq=nnt
 660           nsup=nct-nnt+1
 661           do i=1,2*nres
 662             do j=1,3
 663               cref(j,i)=c(j,i)
 664             enddo
 665           enddo
 666           call contact(.true.,ncont_ref,icont_ref,co)
 667         endif
 668 c        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
 669         call flush(iout)
 670         if (constr_dist.gt.0) call read_dist_constr
 671 c        write (iout,*) "After read_dist_constr nhpb",nhpb
 672         call hpb_partition
 673         if(me.eq.king.or..not.out1file)
 674      &   write (iout,*) 'Contact order:',co
 675         if (pdbref) then
 676         if(me.eq.king.or..not.out1file)
 677      &   write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
 678         do i=1,ncont_ref
 679           do j=1,2
 680             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
 681           enddo
 682           if(me.eq.king.or..not.out1file)
 683      &     write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
 684      &     icont_ref(1,i),' ',
 685      &     restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
 686         enddo
 687         endif
 688       endif
 689       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
 690      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
 691      &    modecalc.ne.10) then
 692 C If input structure hasn't been supplied from the PDB file read or generate
 693 C initial geometry.
 694         if (iranconf.eq.0 .and. .not. extconf) then
 695           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
 696      &     write (iout,'(a)') 'Initial geometry will be read in.'
 697           if (read_cart) then
 698             read(inp,'(8f10.5)',end=36,err=36)
 699      &       ((c(l,k),l=1,3),k=1,nres),
 700      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
 701             call int_from_cart1(.false.)
 702             do i=1,nres-1
 703               do j=1,3
 704                 dc(j,i)=c(j,i+1)-c(j,i)
 705                 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
 706               enddo
 707             enddo
 708             do i=nnt,nct
 709               if (itype(i).ne.10) then
 710                 do j=1,3
 711                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
 712                   dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
 713                 enddo
 714               endif
 715             enddo
 716             return
 717           else
 718             call read_angles(inp,*36)
 719           endif
 720           goto 37
 721    36     write (iout,'(a)') 'Error reading angle file.'
 722 #ifdef MPI
 723           call mpi_finalize( MPI_COMM_WORLD,IERR )
 724 #endif
 725           stop 'Error reading angle file.'
 726    37     continue
 727         else if (extconf) then
 728          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
 729      &    write (iout,'(a)') 'Extended chain initial geometry.'
 730          do i=3,nres
 731           theta(i)=90d0*deg2rad
 732          enddo
 733          do i=4,nres
 734           phi(i)=180d0*deg2rad
 735          enddo
 736          do i=2,nres-1
 737           alph(i)=110d0*deg2rad
 738          enddo
 739          do i=2,nres-1
 740           omeg(i)=-120d0*deg2rad
 741          enddo
 742         else
 743           if(me.eq.king.or..not.out1file)
 744      &     write (iout,'(a)') 'Random-generated initial geometry.'
 745
 746
 747 #ifdef MPI
 748           if (me.eq.king  .or. fg_rank.eq.0 .and. (
 749      &           modecalc.eq.12 .or. modecalc.eq.14) ) then
 750 #endif
 751             do itrial=1,100
 752               itmp=1
 753 c              call gen_rand_conf(itmp,*30)
 754               goto 40
 755    30         write (iout,*) 'Failed to generate random conformation',
 756      &          ', itrial=',itrial
 757               write (*,*) 'Processor:',me,
 758      &          ' Failed to generate random conformation',
 759      &          ' itrial=',itrial
 760               call intout
 761
 762 #ifdef AIX
 763               call flush_(iout)
 764 #else
 765               call flush(iout)
 766 #endif
 767             enddo
 768             write (iout,'(a,i3,a)') 'Processor:',me,
 769      &        ' error in generating random conformation.'
 770             write (*,'(a,i3,a)') 'Processor:',me,
 771      &        ' error in generating random conformation.'
 772             call flush(iout)
 773 #ifdef MPI
 774             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
 775    40       continue
 776           endif
 777 #else
 778           do itrial=1,100
 779             itmp=1
 780 c            call gen_rand_conf(itmp,*31)
 781             goto 40
 782    31       write (iout,*) 'Failed to generate random conformation',
 783      &        ', itrial=',itrial
 784             write (*,*) 'Failed to generate random conformation',
 785      &        ', itrial=',itrial
 786           enddo
 787           write (iout,'(a,i3,a)') 'Processor:',me,
 788      &      ' error in generating random conformation.'
 789           write (*,'(a,i3,a)') 'Processor:',me,
 790      &      ' error in generating random conformation.'
 791           stop
 792    40     continue
 793 #endif
 794         endif
 795       elseif (modecalc.eq.4) then
 796         read (inp,'(a)') intinname
 797         open (intin,file=intinname,status='old',err=333)
 798         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
 799      &  write (iout,'(a)') 'intinname',intinname
 800         write (*,'(a)') 'Processor',myrank,' intinname',intinname
 801         goto 334
 802   333   write (iout,'(2a)') 'Error opening angle file ',intinname
 803 #ifdef MPI
 804         call MPI_Finalize(MPI_COMM_WORLD,IERR)
 805 #endif
 806         stop 'Error opening angle file.'
 807   334   continue
 808
 809       endif
 810 C Generate distance constraints, if the PDB structure is to be regularized.
 811       if (nthread.gt.0) then
 812         call read_threadbase
 813       endif
 814       call setup_var
 815       if (me.eq.king .or. .not. out1file)
 816      & call intout
 817       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
 818         write (iout,'(/a,i3,a)')
 819      &  'The chain contains',ns,' disulfide-bridging cysteines.'
 820         write (iout,'(20i4)') (iss(i),i=1,ns)
 821         write (iout,'(/a/)') 'Pre-formed links are:'
 822         do i=1,nss
 823           i1=ihpb(i)-nres
 824           i2=jhpb(i)-nres
 825           it1=itype(i1)
 826           it2=itype(i2)
 827           if (me.eq.king.or..not.out1file)
 828      &    write (iout,'(2a,i3,3a,i3,a,3f10.3)')
 829      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
 830      &    ebr,forcon(i)
 831         enddo
 832         write (iout,'(a)')
 833       endif
 834 c       if (i2ndstr.gt.0) call secstrp2dihc
 835 c      call geom_to_var(nvar,x)
 836 c      call etotal(energia(0))
 837 c      call enerprint(energia(0))
 838 c      call briefout(0,etot)
 839 c      stop
 840 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
 841 cd    write (iout,'(a)') 'Variable list:'
 842 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
 843 #ifdef MPI
 844       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
 845      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
 846      &  'Processor',myrank,': end reading molecular data.'
 847 #endif
 848       return
 849       end
 850 c--------------------------------------------------------------------------
 851       logical function seq_comp(itypea,itypeb,length)
 852       implicit none
 853       integer length,itypea(length),itypeb(length)
 854       integer i
 855       do i=1,length
 856         if (itypea(i).ne.itypeb(i)) then
 857           seq_comp=.false.
 858           return
 859         endif
 860       enddo
 861       seq_comp=.true.
 862       return
 863       end
 864 c-----------------------------------------------------------------------------
 865       subroutine read_bridge
 866 C Read information about disulfide bridges.
 867       implicit real*8 (a-h,o-z)
 868       include 'DIMENSIONS'
 869 #ifdef MPI
 870       include 'mpif.h'
 871 #endif
 872       include 'COMMON.IOUNITS'
 873       include 'COMMON.GEO'
 874       include 'COMMON.VAR'
 875       include 'COMMON.INTERACT'
 876       include 'COMMON.LOCAL'
 877       include 'COMMON.NAMES'
 878       include 'COMMON.CHAIN'
 879       include 'COMMON.FFIELD'
 880       include 'COMMON.SBRIDGE'
 881       include 'COMMON.HEADER'
 882       include 'COMMON.CONTROL'
 883 c      include 'COMMON.DBASE'
 884 c      include 'COMMON.THREAD'
 885       include 'COMMON.TIME1'
 886       include 'COMMON.SETUP'
 887 C Read bridging residues.
 888       read (inp,*) ns,(iss(i),i=1,ns)
 889       print *,'ns=',ns
 890       if(me.eq.king.or..not.out1file)
 891      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
 892 C Check whether the specified bridging residues are cystines.
 893       do i=1,ns
 894         if (itype(iss(i)).ne.1) then
 895           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
 896      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 897      &   ' can form a disulfide bridge?!!!'
 898           write (*,'(2a,i3,a)')
 899      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 900      &   ' can form a disulfide bridge?!!!'
 901 #ifdef MPI
 902          call MPI_Finalize(MPI_COMM_WORLD,ierror)
 903          stop
 904 #endif
 905         endif
 906       enddo
 907 C Read preformed bridges.
 908       if (ns.gt.0) then
 909       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 910       write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
 911       if (nss.gt.0) then
 912         nhpb=nss
 913 C Check if the residues involved in bridges are in the specified list of
 914 C bridging residues.
 915         do i=1,nss
 916           do j=1,i-1
 917             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 918      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 919               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 920      &      ' contains residues present in other pairs.'
 921               write (*,'(a,i3,a)') 'Disulfide pair',i,
 922      &      ' contains residues present in other pairs.'
 923 #ifdef MPI
 924               call MPI_Finalize(MPI_COMM_WORLD,ierror)
 925               stop
 926 #endif
 927             endif
 928           enddo
 929           do j=1,ns
 930             if (ihpb(i).eq.iss(j)) goto 10
 931           enddo
 932           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 933    10     continue
 934           do j=1,ns
 935             if (jhpb(i).eq.iss(j)) goto 20
 936           enddo
 937           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 938    20     continue
 939           dhpb(i)=dbr
 940           forcon(i)=fbr
 941         enddo
 942         do i=1,nss
 943           ihpb(i)=ihpb(i)+nres
 944           jhpb(i)=jhpb(i)+nres
 945         enddo
 946       endif
 947       endif
 948       return
 949       end
 950 c----------------------------------------------------------------------------
 951       subroutine read_x(kanal,*)
 952       implicit real*8 (a-h,o-z)
 953       include 'DIMENSIONS'
 954       include 'COMMON.GEO'
 955       include 'COMMON.VAR'
 956       include 'COMMON.CHAIN'
 957       include 'COMMON.IOUNITS'
 958       include 'COMMON.CONTROL'
 959       include 'COMMON.LOCAL'
 960       include 'COMMON.INTERACT'
 961 c Read coordinates from input
 962 c
 963       read(kanal,'(8f10.5)',end=10,err=10)
 964      &  ((c(l,k),l=1,3),k=1,nres),
 965      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
 966       do j=1,3
 967         c(j,nres+1)=c(j,1)
 968         c(j,2*nres)=c(j,nres)
 969       enddo
 970       call int_from_cart1(.false.)
 971       do i=1,nres-1
 972         do j=1,3
 973           dc(j,i)=c(j,i+1)-c(j,i)
 974           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 975         enddo
 976       enddo
 977       do i=nnt,nct
 978         if (itype(i).ne.10) then
 979           do j=1,3
 980             dc(j,i+nres)=c(j,i+nres)-c(j,i)
 981             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 982           enddo
 983         endif
 984       enddo
 985
 986       return
 987    10 return1
 988       end
 989 c------------------------------------------------------------------------------
 990       subroutine setup_var
 991       implicit real*8 (a-h,o-z)
 992       include 'DIMENSIONS'
 993       include 'COMMON.IOUNITS'
 994       include 'COMMON.GEO'
 995       include 'COMMON.VAR'
 996       include 'COMMON.INTERACT'
 997       include 'COMMON.LOCAL'
 998       include 'COMMON.NAMES'
 999       include 'COMMON.CHAIN'
1000       include 'COMMON.FFIELD'
1001       include 'COMMON.SBRIDGE'
1002       include 'COMMON.HEADER'
1003       include 'COMMON.CONTROL'
1004 c      include 'COMMON.DBASE'
1005 c      include 'COMMON.THREAD'
1006       include 'COMMON.TIME1'
1007 C Set up variable list.
1008       ntheta=nres-2
1009       nphi=nres-3
1010       nvar=ntheta+nphi
1011       nside=0
1012       do i=2,nres-1
1013         if (itype(i).ne.10) then
1014           nside=nside+1
1015           ialph(i,1)=nvar+nside
1016           ialph(nside,2)=i
1017         endif
1018       enddo
1019       if (indphi.gt.0) then
1020         nvar=nphi
1021       else if (indback.gt.0) then
1022         nvar=nphi+ntheta
1023       else
1024         nvar=nvar+2*nside
1025       endif
1026 cd    write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1027       return
1028       end
1029 c----------------------------------------------------------------------------
1030       subroutine gen_dist_constr
1031 C Generate CA distance constraints.
1032       implicit real*8 (a-h,o-z)
1033       include 'DIMENSIONS'
1034       include 'COMMON.IOUNITS'
1035       include 'COMMON.GEO'
1036       include 'COMMON.VAR'
1037       include 'COMMON.INTERACT'
1038       include 'COMMON.LOCAL'
1039       include 'COMMON.NAMES'
1040       include 'COMMON.CHAIN'
1041       include 'COMMON.FFIELD'
1042       include 'COMMON.SBRIDGE'
1043       include 'COMMON.HEADER'
1044       include 'COMMON.CONTROL'
1045 c      include 'COMMON.DBASE'
1046 c      include 'COMMON.THREAD'
1047       include 'COMMON.TIME1'
1048       dimension itype_pdb(maxres)
1049       common /pizda/ itype_pdb
1050       character*2 iden
1051 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1052 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1053 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1054 cd     & ' nsup',nsup
1055       do i=nstart_sup,nstart_sup+nsup-1
1056 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1057 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1058         do j=i+2,nstart_sup+nsup-1
1059           nhpb=nhpb+1
1060           ihpb(nhpb)=i+nstart_seq-nstart_sup
1061           jhpb(nhpb)=j+nstart_seq-nstart_sup
1062           forcon(nhpb)=weidis
1063           dhpb(nhpb)=dist(i,j)
1064         enddo
1065       enddo
1066 cd      write (iout,'(a)') 'Distance constraints:'
1067 cd      do i=nss+1,nhpb
1068 cd        ii=ihpb(i)
1069 cd        jj=jhpb(i)
1070 cd        iden='CA'
1071 cd        if (ii.gt.nres) then
1072 cd          iden='SC'
1073 cd          ii=ii-nres
1074 cd          jj=jj-nres
1075 cd        endif
1076 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1077 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1078 cd     &  dhpb(i),forcon(i)
1079 cd      enddo
1080       return
1081       end
1082 c----------------------------------------------------------------------------
1083       subroutine csaread
1084       implicit real*8 (a-h,o-z)
1085       include 'DIMENSIONS'
1086       include 'COMMON.IOUNITS'
1087       include 'COMMON.GEO'
1088       include 'COMMON.CSA'
1089       include 'COMMON.BANK'
1090       include 'COMMON.CONTROL'
1091       character*80 ucase
1092       character*620 mcmcard
1093       call card_concat(mcmcard)
1094
1095       call readi(mcmcard,'NCONF',nconf,50)
1096       call readi(mcmcard,'NADD',nadd,0)
1097       call readi(mcmcard,'JSTART',jstart,1)
1098       call readi(mcmcard,'JEND',jend,1)
1099       call readi(mcmcard,'NSTMAX',nstmax,500000)
1100       call readi(mcmcard,'N0',n0,1)
1101       call readi(mcmcard,'N1',n1,6)
1102       call readi(mcmcard,'N2',n2,4)
1103       call readi(mcmcard,'N3',n3,0)
1104       call readi(mcmcard,'N4',n4,0)
1105       call readi(mcmcard,'N5',n5,0)
1106       call readi(mcmcard,'N6',n6,10)
1107       call readi(mcmcard,'N7',n7,0)
1108       call readi(mcmcard,'N8',n8,0)
1109       call readi(mcmcard,'N9',n9,0)
1110       call readi(mcmcard,'N14',n14,0)
1111       call readi(mcmcard,'N15',n15,0)
1112       call readi(mcmcard,'N16',n16,0)
1113       call readi(mcmcard,'N17',n17,0)
1114       call readi(mcmcard,'N18',n18,0)
1115
1116       vdisulf=(index(mcmcard,'DYNSS').gt.0)
1117
1118       call readi(mcmcard,'NDIFF',ndiff,2)
1119       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1120       call readi(mcmcard,'IS1',is1,1)
1121       call readi(mcmcard,'IS2',is2,8)
1122       call readi(mcmcard,'NRAN0',nran0,4)
1123       call readi(mcmcard,'NRAN1',nran1,2)
1124       call readi(mcmcard,'IRR',irr,1)
1125       call readi(mcmcard,'NSEED',nseed,20)
1126       call readi(mcmcard,'NTOTAL',ntotal,10000)
1127       call reada(mcmcard,'CUT1',cut1,2.0d0)
1128       call reada(mcmcard,'CUT2',cut2,5.0d0)
1129       call reada(mcmcard,'ESTOP',estop,-300000.0d0)
1130       call readi(mcmcard,'ICMAX',icmax,1)
1131       call readi(mcmcard,'NBANKM',nbankm,400)
1132       call readi(mcmcard,'IUCUT',iucut,2)
1133       call readi(mcmcard,'IRESTART',irestart,0)
1134 c!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
1135       ntbankm=0
1136 c!bankt
1137       call reada(mcmcard,'DELE',dele,20.0d0)
1138       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1139       call readi(mcmcard,'IREF',iref,0)
1140       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1141       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1142       call readi(mcmcard,'NCONF_IN',nconf_in,0)
1143       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1144       write (iout,*) "NCONF_IN",nconf_in
1145       tm_score=(index(mcmcard,'TMSCORE').gt.0)
1146       if (tm_score) write (iout,*) "Using TM_Score instead of DIFF",
1147      &      " for torsional angles"
1148       return
1149       end
1150
1151       subroutine read_minim
1152       implicit real*8 (a-h,o-z)
1153       include 'DIMENSIONS'
1154       include 'COMMON.MINIM'
1155       include 'COMMON.IOUNITS'
1156       character*80 ucase
1157       character*320 minimcard
1158       call card_concat(minimcard)
1159       call readi(minimcard,'MAXMIN',maxmin,2000)
1160       call readi(minimcard,'MAXFUN',maxfun,5000)
1161       call readi(minimcard,'MINMIN',minmin,maxmin)
1162       call readi(minimcard,'MINFUN',minfun,maxmin)
1163       call reada(minimcard,'TOLF',tolf,1.0D-2)
1164       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1165       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1166       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1167       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1168       write (iout,'(/80(1h*)/20x,a/80(1h*))')
1169      &         'Options in energy minimization:'
1170       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1171      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1172      & 'MinMin:',MinMin,' MinFun:',MinFun,
1173      & ' TolF:',TolF,' RTolF:',RTolF
1174       return
1175       end
1176 c----------------------------------------------------------------------------
1177       subroutine read_angles(kanal,*)
1178       implicit real*8 (a-h,o-z)
1179       include 'DIMENSIONS'
1180       include 'COMMON.GEO'
1181       include 'COMMON.VAR'
1182       include 'COMMON.CHAIN'
1183       include 'COMMON.IOUNITS'
1184       include 'COMMON.CONTROL'
1185 c Read angles from input
1186 c
1187        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1188        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1189        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1190        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1191
1192        do i=1,nres
1193 c 9/7/01 avoid 180 deg valence angle
1194         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1195 c
1196         theta(i)=deg2rad*theta(i)
1197         phi(i)=deg2rad*phi(i)
1198         alph(i)=deg2rad*alph(i)
1199         omeg(i)=deg2rad*omeg(i)
1200        enddo
1201       return
1202    10 return1
1203       end
1204 c----------------------------------------------------------------------------
1205       subroutine reada(rekord,lancuch,wartosc,default)
1206       implicit none
1207       character*(*) rekord,lancuch
1208       double precision wartosc,default
1209       integer ilen,iread
1210       external ilen
1211       iread=index(rekord,lancuch)
1212       if (iread.eq.0) then
1213         wartosc=default
1214         return
1215       endif
1216       iread=iread+ilen(lancuch)+1
1217       read (rekord(iread:),*,err=10,end=10) wartosc
1218       return
1219   10  wartosc=default
1220       return
1221       end
1222 c----------------------------------------------------------------------------
1223       subroutine readi(rekord,lancuch,wartosc,default)
1224       implicit none
1225       character*(*) rekord,lancuch
1226       integer wartosc,default
1227       integer ilen,iread
1228       external ilen
1229       iread=index(rekord,lancuch)
1230       if (iread.eq.0) then
1231         wartosc=default
1232         return
1233       endif
1234       iread=iread+ilen(lancuch)+1
1235       read (rekord(iread:),*,err=10,end=10) wartosc
1236       return
1237   10  wartosc=default
1238       return
1239       end
1240 c----------------------------------------------------------------------------
1241       subroutine multreadi(rekord,lancuch,tablica,dim,default)
1242       implicit none
1243       integer dim,i
1244       integer tablica(dim),default
1245       character*(*) rekord,lancuch
1246       character*80 aux
1247       integer ilen,iread
1248       external ilen
1249       do i=1,dim
1250         tablica(i)=default
1251       enddo
1252       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1253       if (iread.eq.0) return
1254       iread=iread+ilen(lancuch)+1
1255       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1256    10 return
1257       end
1258 c----------------------------------------------------------------------------
1259       subroutine multreada(rekord,lancuch,tablica,dim,default)
1260       implicit none
1261       integer dim,i
1262       double precision tablica(dim),default
1263       character*(*) rekord,lancuch
1264       character*80 aux
1265       integer ilen,iread
1266       external ilen
1267       do i=1,dim
1268         tablica(i)=default
1269       enddo
1270       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1271       if (iread.eq.0) return
1272       iread=iread+ilen(lancuch)+1
1273       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1274    10 return
1275       end
1276 c----------------------------------------------------------------------------
1277       subroutine openunits
1278       implicit real*8 (a-h,o-z)
1279       include 'DIMENSIONS'
1280 #ifdef MPI
1281       include 'mpif.h'
1282       character*16 form,nodename
1283       integer nodelen
1284 #endif
1285       include 'COMMON.SETUP'
1286       include 'COMMON.IOUNITS'
1287 c      include 'COMMON.MD'
1288       include 'COMMON.CONTROL'
1289       integer lenpre,lenpot,ilen,lentmp
1290       external ilen
1291       character*3 out1file_text,ucase
1292       character*3 ll
1293       external ucase
1294 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1295       call getenv_loc("PREFIX",prefix)
1296       pref_orig = prefix
1297       call getenv_loc("POT",pot)
1298       call getenv_loc("DIRTMP",tmpdir)
1299       call getenv_loc("CURDIR",curdir)
1300       call getenv_loc("OUT1FILE",out1file_text)
1301 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1302       out1file_text=ucase(out1file_text)
1303       if (out1file_text(1:1).eq."Y") then
1304         out1file=.true.
1305       else
1306         out1file=fg_rank.gt.0
1307       endif
1308       lenpre=ilen(prefix)
1309       lenpot=ilen(pot)
1310       lentmp=ilen(tmpdir)
1311       if (lentmp.gt.0) then
1312           write (*,'(80(1h!))')
1313           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
1314           write (*,'(80(1h!))')
1315           write (*,*)"All output files will be on node /tmp directory."
1316 #ifdef MPI
1317         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1318         if (me.eq.king) then
1319           write (*,*) "The master node is ",nodename
1320         else if (fg_rank.eq.0) then
1321           write (*,*) "I am the CG slave node ",nodename
1322         else
1323           write (*,*) "I am the FG slave node ",nodename
1324         endif
1325 #endif
1326         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1327         lenpre = lentmp+lenpre+1
1328       endif
1329       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1330 C Get the names and open the input files
1331 #if defined(WINIFL) || defined(WINPGI)
1332       open(1,file=pref_orig(:ilen(pref_orig))//
1333      &  '.inp',status='old',readonly,shared)
1334        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1335 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1336 C Get parameter filenames and open the parameter files.
1337       call getenv_loc('BONDPAR',bondname)
1338       open (ibond,file=bondname,status='old',readonly,shared)
1339       call getenv_loc('THETPAR',thetname)
1340       open (ithep,file=thetname,status='old',readonly,shared)
1341 #ifndef CRYST_THETA
1342       call getenv_loc('THETPARPDB',thetname_pdb)
1343       open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
1344 #endif
1345       call getenv_loc('ROTPAR',rotname)
1346       open (irotam,file=rotname,status='old',readonly,shared)
1347 #ifndef CRYST_SC
1348       call getenv_loc('ROTPARPDB',rotname_pdb)
1349       open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
1350 #endif
1351       call getenv_loc('TORPAR',torname)
1352       open (itorp,file=torname,status='old',readonly,shared)
1353       call getenv_loc('TORDPAR',tordname)
1354       open (itordp,file=tordname,status='old',readonly,shared)
1355       call getenv_loc('FOURIER',fouriername)
1356       open (ifourier,file=fouriername,status='old',readonly,shared)
1357       call getenv_loc('ELEPAR',elename)
1358       open (ielep,file=elename,status='old',readonly,shared)
1359       call getenv_loc('SIDEPAR',sidename)
1360       open (isidep,file=sidename,status='old',readonly,shared)
1361 #elif (defined CRAY) || (defined AIX)
1362       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1363      &  action='read')
1364 c      print *,"Processor",myrank," opened file 1"
1365       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1366 c      print *,"Processor",myrank," opened file 9"
1367 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1368 C Get parameter filenames and open the parameter files.
1369       call getenv_loc('BONDPAR',bondname)
1370       open (ibond,file=bondname,status='old',action='read')
1371 c      print *,"Processor",myrank," opened file IBOND"
1372       call getenv_loc('THETPAR',thetname)
1373       open (ithep,file=thetname,status='old',action='read')
1374 c      print *,"Processor",myrank," opened file ITHEP"
1375 #ifndef CRYST_THETA
1376       call getenv_loc('THETPARPDB',thetname_pdb)
1377       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1378 #endif
1379       call getenv_loc('ROTPAR',rotname)
1380       open (irotam,file=rotname,status='old',action='read')
1381 c      print *,"Processor",myrank," opened file IROTAM"
1382 #ifndef CRYST_SC
1383       call getenv_loc('ROTPARPDB',rotname_pdb)
1384       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
1385 #endif
1386       call getenv_loc('TORPAR',torname)
1387       open (itorp,file=torname,status='old',action='read')
1388 c      print *,"Processor",myrank," opened file ITORP"
1389       call getenv_loc('TORDPAR',tordname)
1390       open (itordp,file=tordname,status='old',action='read')
1391 c      print *,"Processor",myrank," opened file ITORDP"
1392       call getenv_loc('SCCORPAR',sccorname)
1393       open (isccor,file=sccorname,status='old',action='read')
1394 c      print *,"Processor",myrank," opened file ISCCOR"
1395       call getenv_loc('FOURIER',fouriername)
1396       open (ifourier,file=fouriername,status='old',action='read')
1397 c      print *,"Processor",myrank," opened file IFOURIER"
1398       call getenv_loc('ELEPAR',elename)
1399       open (ielep,file=elename,status='old',action='read')
1400 c      print *,"Processor",myrank," opened file IELEP"
1401       call getenv_loc('SIDEPAR',sidename)
1402       open (isidep,file=sidename,status='old',action='read')
1403 c      print *,"Processor",myrank," opened file ISIDEP"
1404 c      print *,"Processor",myrank," opened parameter files"
1405 #elif (defined G77)
1406       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1407       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1408 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1409 C Get parameter filenames and open the parameter files.
1410       call getenv_loc('BONDPAR',bondname)
1411       open (ibond,file=bondname,status='old')
1412       call getenv_loc('THETPAR',thetname)
1413       open (ithep,file=thetname,status='old')
1414 #ifndef CRYST_THETA
1415       call getenv_loc('THETPARPDB',thetname_pdb)
1416       open (ithep_pdb,file=thetname_pdb,status='old')
1417 #endif
1418       call getenv_loc('ROTPAR',rotname)
1419       open (irotam,file=rotname,status='old')
1420 #ifndef CRYST_SC
1421       call getenv_loc('ROTPARPDB',rotname_pdb)
1422       open (irotam_pdb,file=rotname_pdb,status='old')
1423 #endif
1424       call getenv_loc('TORPAR',torname)
1425       open (itorp,file=torname,status='old')
1426       call getenv_loc('TORDPAR',tordname)
1427       open (itordp,file=tordname,status='old')
1428       call getenv_loc('SCCORPAR',sccorname)
1429       open (isccor,file=sccorname,status='old')
1430       call getenv_loc('FOURIER',fouriername)
1431       open (ifourier,file=fouriername,status='old')
1432       call getenv_loc('ELEPAR',elename)
1433       open (ielep,file=elename,status='old')
1434       call getenv_loc('SIDEPAR',sidename)
1435       open (isidep,file=sidename,status='old')
1436 #else
1437       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1438      &  readonly)
1439        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1440 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1441 C Get parameter filenames and open the parameter files.
1442       call getenv_loc('BONDPAR',bondname)
1443       open (ibond,file=bondname,status='old',readonly)
1444       call getenv_loc('THETPAR',thetname)
1445       open (ithep,file=thetname,status='old',readonly)
1446 #ifndef CRYST_THETA
1447       call getenv_loc('THETPARPDB',thetname_pdb)
1448       print *,"thetname_pdb ",thetname_pdb
1449       open (ithep_pdb,file=thetname_pdb,status='old',readonly)
1450       print *,ithep_pdb," opened"
1451 #endif
1452       call getenv_loc('ROTPAR',rotname)
1453       open (irotam,file=rotname,status='old',readonly)
1454 #ifndef CRYST_SC
1455       call getenv_loc('ROTPARPDB',rotname_pdb)
1456       open (irotam_pdb,file=rotname_pdb,status='old',readonly)
1457 #endif
1458       call getenv_loc('TORPAR',torname)
1459       open (itorp,file=torname,status='old',readonly)
1460       call getenv_loc('TORDPAR',tordname)
1461       open (itordp,file=tordname,status='old',readonly)
1462       call getenv_loc('SCCORPAR',sccorname)
1463       open (isccor,file=sccorname,status='old',readonly)
1464       call getenv_loc('FOURIER',fouriername)
1465       open (ifourier,file=fouriername,status='old',readonly)
1466       call getenv_loc('ELEPAR',elename)
1467       open (ielep,file=elename,status='old',readonly)
1468       call getenv_loc('SIDEPAR',sidename)
1469       open (isidep,file=sidename,status='old',readonly)
1470 #endif
1471 #ifndef OLDSCP
1472 C
1473 C 8/9/01 In the newest version SCp interaction constants are read from a file
1474 C Use -DOLDSCP to use hard-coded constants instead.
1475 C
1476       call getenv_loc('SCPPAR',scpname)
1477 #if defined(WINIFL) || defined(WINPGI)
1478       open (iscpp,file=scpname,status='old',readonly,shared)
1479 #elif (defined CRAY)  || (defined AIX)
1480       open (iscpp,file=scpname,status='old',action='read')
1481 #elif (defined G77)
1482       open (iscpp,file=scpname,status='old')
1483 #else
1484       open (iscpp,file=scpname,status='old',readonly)
1485 #endif
1486 #endif
1487       call getenv_loc('PATTERN',patname)
1488 #if defined(WINIFL) || defined(WINPGI)
1489       open (icbase,file=patname,status='old',readonly,shared)
1490 #elif (defined CRAY)  || (defined AIX)
1491       open (icbase,file=patname,status='old',action='read')
1492 #elif (defined G77)
1493       open (icbase,file=patname,status='old')
1494 #else
1495       open (icbase,file=patname,status='old',readonly)
1496 #endif
1497 #ifdef MPI
1498 C Open output file only for CG processes
1499 c      print *,"Processor",myrank," fg_rank",fg_rank
1500       if (fg_rank.eq.0) then
1501
1502       if (nodes.eq.1) then
1503         npos=3
1504       else
1505         npos = dlog10(dfloat(nodes-1))+1
1506       endif
1507       if (npos.lt.3) npos=3
1508       write (liczba,'(i1)') npos
1509       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
1510      &  //')'
1511       write (liczba,form) me
1512       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
1513      &  liczba(:ilen(liczba))
1514       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1515      &  //'.int'
1516       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1517      &  //'.pdb'
1518       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
1519      &  liczba(:ilen(liczba))//'.mol2'
1520       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
1521      &  liczba(:ilen(liczba))//'.stat'
1522       if (lentmp.gt.0)
1523      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
1524      &      //liczba(:ilen(liczba))//'.stat')
1525       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
1526      &  //'.rst'
1527 c      if(usampl) then
1528 c          qname=prefix(:lenpre)//'_'//pot(:lenpot)//
1529 c     & liczba(:ilen(liczba))//'.const'
1530 c      endif
1531
1532       endif
1533 #else
1534       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
1535       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
1536       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
1537       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
1538       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
1539       if (lentmp.gt.0)
1540      & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)//
1541      &    '.stat')
1542       rest2name=prefix(:ilen(prefix))//'.rst'
1543 c      if(usampl) then
1544 c         qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
1545 c      endif
1546 #endif
1547 #if defined(AIX) || defined(PGI)
1548       if (me.eq.king .or. .not. out1file)
1549      &   open(iout,file=outname,status='unknown')
1550 #ifdef DEBUG
1551       if (fg_rank.gt.0) then
1552         write (liczba,'(i3.3)') myrank/nfgtasks
1553         write (ll,'(bz,i3.3)') fg_rank
1554         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1555      &   status='unknown')
1556       endif
1557 #endif
1558       if(me.eq.king) then
1559        open(igeom,file=intname,status='unknown',position='append')
1560        open(ipdb,file=pdbname,status='unknown')
1561        open(imol2,file=mol2name,status='unknown')
1562        open(istat,file=statname,status='unknown',position='append')
1563       else
1564 c1out       open(iout,file=outname,status='unknown')
1565       endif
1566 #else
1567       if (me.eq.king .or. .not.out1file)
1568      &    open(iout,file=outname,status='unknown')
1569 #ifdef DEBUG
1570       if (fg_rank.gt.0) then
1571         write (liczba,'(i3.3)') myrank/nfgtasks
1572         write (ll,'(bz,i3.3)') fg_rank
1573         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1574      &   status='unknown')
1575       endif
1576 #endif
1577       if(me.eq.king) then
1578        open(igeom,file=intname,status='unknown',access='append')
1579        open(ipdb,file=pdbname,status='unknown')
1580        open(imol2,file=mol2name,status='unknown')
1581        open(istat,file=statname,status='unknown',access='append')
1582       else
1583 c1out       open(iout,file=outname,status='unknown')
1584       endif
1585 #endif
1586       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
1587       csa_seed=prefix(:lenpre)//'.CSA.seed'
1588       csa_history=prefix(:lenpre)//'.CSA.history'
1589       csa_bank=prefix(:lenpre)//'.CSA.bank'
1590       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
1591       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
1592       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
1593 c!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
1594       csa_int=prefix(:lenpre)//'.int'
1595       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
1596       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
1597       csa_in=prefix(:lenpre)//'.CSA.in'
1598 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
1599 C Write file names
1600       if (me.eq.king)then
1601       write (iout,'(80(1h-))')
1602       write (iout,'(30x,a)') "FILE ASSIGNMENT"
1603       write (iout,'(80(1h-))')
1604       write (iout,*) "Input file                      : ",
1605      &  pref_orig(:ilen(pref_orig))//'.inp'
1606       write (iout,*) "Output file                     : ",
1607      &  outname(:ilen(outname))
1608       write (iout,*)
1609       write (iout,*) "Sidechain potential file        : ",
1610      &  sidename(:ilen(sidename))
1611 #ifndef OLDSCP
1612       write (iout,*) "SCp potential file              : ",
1613      &  scpname(:ilen(scpname))
1614 #endif
1615       write (iout,*) "Electrostatic potential file    : ",
1616      &  elename(:ilen(elename))
1617       write (iout,*) "Cumulant coefficient file       : ",
1618      &  fouriername(:ilen(fouriername))
1619       write (iout,*) "Torsional parameter file        : ",
1620      &  torname(:ilen(torname))
1621       write (iout,*) "Double torsional parameter file : ",
1622      &  tordname(:ilen(tordname))
1623       write (iout,*) "SCCOR parameter file : ",
1624      &  sccorname(:ilen(sccorname))
1625       write (iout,*) "Bond & inertia constant file    : ",
1626      &  bondname(:ilen(bondname))
1627       write (iout,*) "Bending parameter file          : ",
1628      &  thetname(:ilen(thetname))
1629       write (iout,*) "Rotamer parameter file          : ",
1630      &  rotname(:ilen(rotname))
1631       write (iout,*) "Threading database              : ",
1632      &  patname(:ilen(patname))
1633       if (lentmp.ne.0)
1634      &write (iout,*)" DIRTMP                          : ",
1635      &  tmpdir(:lentmp)
1636       write (iout,'(80(1h-))')
1637       endif
1638       return
1639       end
1640 c----------------------------------------------------------------------------
1641       subroutine card_concat(card)
1642       implicit real*8 (a-h,o-z)
1643       include 'DIMENSIONS'
1644       include 'COMMON.IOUNITS'
1645       character*(*) card
1646       character*80 karta,ucase
1647       external ilen
1648       read (inp,'(a)') karta
1649       karta=ucase(karta)
1650       card=' '
1651       do while (karta(80:80).eq.'&')
1652         card=card(:ilen(card)+1)//karta(:79)
1653         read (inp,'(a)') karta
1654         karta=ucase(karta)
1655       enddo
1656       card=card(:ilen(card)+1)//karta
1657       return
1658       end
1659
1660       subroutine read_dist_constr
1661       implicit real*8 (a-h,o-z)
1662       include 'DIMENSIONS'
1663 #ifdef MPI
1664       include 'mpif.h'
1665 #endif
1666       include 'COMMON.SETUP'
1667       include 'COMMON.CONTROL'
1668       include 'COMMON.CHAIN'
1669       include 'COMMON.IOUNITS'
1670       include 'COMMON.SBRIDGE'
1671       integer ifrag_(2,100),ipair_(2,100)
1672       double precision wfrag_(100),wpair_(100)
1673       character*500 controlcard
1674       write (iout,*) "Calling read_dist_constr"
1675       write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
1676       call flush(iout)
1677       call card_concat(controlcard)
1678       call readi(controlcard,"NFRAG",nfrag_,0)
1679       call readi(controlcard,"NPAIR",npair_,0)
1680       call readi(controlcard,"NDIST",ndist_,0)
1681       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
1682       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
1683       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
1684       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
1685       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
1686 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
1687 c      write (iout,*) "IFRAG"
1688 c      do i=1,nfrag_
1689 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1690 c      enddo
1691 c      write (iout,*) "IPAIR"
1692 c      do i=1,npair_
1693 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
1694 c      enddo
1695       call flush(iout)
1696       do i=1,nfrag_
1697         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
1698         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
1699      &    ifrag_(2,i)=nstart_sup+nsup-1
1700 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1701         call flush(iout)
1702         if (wfrag_(i).gt.0.0d0) then
1703         do j=ifrag_(1,i),ifrag_(2,i)-1
1704           do k=j+1,ifrag_(2,i)
1705             write (iout,*) "j",j," k",k
1706             ddjk=dist(j,k)
1707             if (constr_dist.eq.1) then
1708             nhpb=nhpb+1
1709             ihpb(nhpb)=j
1710             jhpb(nhpb)=k
1711               dhpb(nhpb)=ddjk
1712             forcon(nhpb)=wfrag_(i)
1713             else if (constr_dist.eq.2) then
1714               if (ddjk.le.dist_cut) then
1715                 nhpb=nhpb+1
1716                 ihpb(nhpb)=j
1717                 jhpb(nhpb)=k
1718                 dhpb(nhpb)=ddjk
1719                 forcon(nhpb)=wfrag_(i)
1720               endif
1721             else
1722               nhpb=nhpb+1
1723               ihpb(nhpb)=j
1724               jhpb(nhpb)=k
1725               dhpb(nhpb)=ddjk
1726               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1727             endif
1728 #ifdef MPI
1729             if (.not.out1file .or. me.eq.king)
1730      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1731      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1732 #else
1733             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1734      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1735 #endif
1736           enddo
1737         enddo
1738         endif
1739       enddo
1740       do i=1,npair_
1741         if (wpair_(i).gt.0.0d0) then
1742         ii = ipair_(1,i)
1743         jj = ipair_(2,i)
1744         if (ii.gt.jj) then
1745           itemp=ii
1746           ii=jj
1747           jj=itemp
1748         endif
1749         do j=ifrag_(1,ii),ifrag_(2,ii)
1750           do k=ifrag_(1,jj),ifrag_(2,jj)
1751             nhpb=nhpb+1
1752             ihpb(nhpb)=j
1753             jhpb(nhpb)=k
1754             forcon(nhpb)=wpair_(i)
1755             dhpb(nhpb)=dist(j,k)
1756 #ifdef MPI
1757             if (.not.out1file .or. me.eq.king)
1758      &      write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1759      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1760 #else
1761             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1762      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1763 #endif
1764           enddo
1765         enddo
1766         endif
1767       enddo
1768       do i=1,ndist_
1769         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
1770         if (forcon(nhpb+1).gt.0.0d0) then
1771           nhpb=nhpb+1
1772           dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1773 #ifdef MPI
1774           if (.not.out1file .or. me.eq.king)
1775      &    write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1776      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1777 #else
1778           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1779      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1780 #endif
1781         endif
1782       enddo
1783       call flush(iout)
1784       return
1785       end
1786 c-------------------------------------------------------------------------------
1787 #ifdef WINIFL
1788       subroutine flush(iu)
1789       return
1790       end
1791 #endif
1792 #ifdef AIX
1793       subroutine flush(iu)
1794       call flush_(iu)
1795       return
1796       end
1797 #endif
1798 c------------------------------------------------------------------------------
1799       subroutine copy_to_tmp(source)
1800       include "DIMENSIONS"
1801       include "COMMON.IOUNITS"
1802       character*(*) source
1803       character* 256 tmpfile
1804       integer ilen
1805       external ilen
1806       logical ex
1807       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
1808       inquire(file=tmpfile,exist=ex)
1809       if (ex) then
1810         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
1811      &   " to temporary directory..."
1812         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
1813         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
1814       endif
1815       return
1816       end
1817 c------------------------------------------------------------------------------
1818       subroutine move_from_tmp(source)
1819       include "DIMENSIONS"
1820       include "COMMON.IOUNITS"
1821       character*(*) source
1822       integer ilen
1823       external ilen
1824       write (*,*) "Moving ",source(:ilen(source)),
1825      & " from temporary directory to working directory"
1826       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
1827       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
1828       return
1829       end
1830 c------------------------------------------------------------------------------
1831       subroutine random_init(seed)
1832 C
1833 C Initialize random number generator
1834 C
1835       implicit real*8 (a-h,o-z)
1836       include 'DIMENSIONS'
1837 #ifdef AMD64
1838       integer*8 iseedi8
1839 #endif
1840 #ifdef MPI
1841       include 'mpif.h'
1842       logical OKRandom, prng_restart
1843       real*8  r1
1844       integer iseed_array(4)
1845 #endif
1846       include 'COMMON.IOUNITS'
1847       include 'COMMON.TIME1'
1848 c      include 'COMMON.THREAD'
1849       include 'COMMON.SBRIDGE'
1850       include 'COMMON.CONTROL'
1851       include 'COMMON.MCM'
1852 c      include 'COMMON.MAP'
1853       include 'COMMON.HEADER'
1854 c      include 'COMMON.CSA'
1855       include 'COMMON.CHAIN'
1856 c      include 'COMMON.MUCA'
1857 c      include 'COMMON.MD'
1858       include 'COMMON.FFIELD'
1859       include 'COMMON.SETUP'
1860       iseed=-dint(dabs(seed))
1861       if (iseed.eq.0) then
1862         write (iout,'(/80(1h*)/20x,a/80(1h*))')
1863      &    'Random seed undefined. The program will stop.'
1864         write (*,'(/80(1h*)/20x,a/80(1h*))')
1865      &    'Random seed undefined. The program will stop.'
1866 #ifdef MPI
1867         call mpi_finalize(mpi_comm_world,ierr)
1868 #endif
1869         stop 'Bad random seed.'
1870       endif
1871 #ifdef MPI
1872       if (fg_rank.eq.0) then
1873       seed=seed*(me+1)+1
1874 #ifdef AMD64
1875       iseedi8=dint(seed)
1876       if(me.eq.king .or. .not. out1file)
1877      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1878       write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1879       OKRandom = prng_restart(me,iseedi8)
1880 #else
1881       do i=1,4
1882        tmp=65536.0d0**(4-i)
1883        iseed_array(i) = dint(seed/tmp)
1884        seed=seed-iseed_array(i)*tmp
1885       enddo
1886       if(me.eq.king .or. .not. out1file)
1887      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
1888      &                 (iseed_array(i),i=1,4)
1889       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
1890      &                 (iseed_array(i),i=1,4)
1891       OKRandom = prng_restart(me,iseed_array)
1892 #endif
1893       if (OKRandom) then
1894         r1=ran_number(0.0D0,1.0D0)
1895         if(me.eq.king .or. .not. out1file)
1896      &   write (iout,*) 'ran_num',r1
1897         if (r1.lt.0.0d0) OKRandom=.false.
1898       endif
1899       if (.not.OKRandom) then
1900         write (iout,*) 'PRNG IS NOT WORKING!!!'
1901         print *,'PRNG IS NOT WORKING!!!'
1902         if (me.eq.0) then
1903          call flush(iout)
1904          call mpi_abort(mpi_comm_world,error_msg,ierr)
1905          stop
1906         else
1907          write (iout,*) 'too many processors for parallel prng'
1908          write (*,*) 'too many processors for parallel prng'
1909          call flush(iout)
1910          stop
1911         endif
1912       endif
1913       endif
1914 #else
1915       call vrndst(iseed)
1916       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
1917 #endif
1918       return
1919       end