2 implicit real*8 (a-h,o-z)
8 include 'COMMON.CONTROL'
9 include 'COMMON.SBRIDGE'
10 include 'COMMON.IOUNITS'
12 C Read force-field parameters except weights
14 C Read job setup parameters
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 if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
20 C Read MD control parameters if reqjuired
21 if (modecalc.eq.12) call read_MDpar
22 C Read MREMD control parameters if required
23 if (modecalc.eq.14) then
27 C Read MUCA control parameters if required
28 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
35 cfmc if (modecalc.eq.10) call mcmfread
36 C Read molecule information, molecule geometry, energy-term weights, and
37 C restraints if requested
39 C Print restraint information
41 if (.not. out1file .or. me.eq.king) then
44 &write (iout,'(a,i5,a)') "The following",nhpb-nss,
45 & " distance constraints have been imposed"
47 write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
52 c print *,"Processor",myrank," leaves READRTNS"
55 C-------------------------------------------------------------------------------
56 subroutine read_control
60 implicit real*8 (a-h,o-z)
64 logical OKRandom, prng_restart
67 include 'COMMON.IOUNITS'
68 include 'COMMON.TIME1'
69 include 'COMMON.THREAD'
70 include 'COMMON.SBRIDGE'
71 include 'COMMON.CONTROL'
74 include 'COMMON.HEADER'
76 include 'COMMON.CHAIN'
79 include 'COMMON.FFIELD'
80 include 'COMMON.INTERACT'
81 include 'COMMON.SETUP'
82 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
83 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
85 character*320 controlcard
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 C this variable with_theta_constr is the variable which allow to read and execute the
100 C constrains on theta angles WITH_THETA_CONSTR is the keyword
101 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
102 write (iout,*) "with_theta_constr ",with_theta_constr
103 call readi(controlcard,'SYM',symetr,1)
104 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
105 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
106 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
107 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
108 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
109 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
110 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
111 call reada(controlcard,'DRMS',drms,0.1D0)
112 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
113 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
114 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
115 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
116 write (iout,'(a,f10.1)')'DRMS = ',drms
117 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
118 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
120 call readi(controlcard,'NZ_START',nz_start,0)
121 call readi(controlcard,'NZ_END',nz_end,0)
122 call readi(controlcard,'IZ_SC',iz_sc,0)
124 safety = 60.0d0*safety
127 call reada(controlcard,"T_BATH",t_bath,300.0d0)
128 minim=(index(controlcard,'MINIMIZE').gt.0)
129 dccart=(index(controlcard,'CART').gt.0)
130 overlapsc=(index(controlcard,'OVERLAP').gt.0)
131 overlapsc=.not.overlapsc
132 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
133 searchsc=.not.searchsc
134 sideadd=(index(controlcard,'SIDEADD').gt.0)
135 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
136 outpdb=(index(controlcard,'PDBOUT').gt.0)
137 outmol2=(index(controlcard,'MOL2OUT').gt.0)
138 pdbref=(index(controlcard,'PDBREF').gt.0)
139 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
140 indpdb=index(controlcard,'PDBSTART')
141 extconf=(index(controlcard,'EXTCONF').gt.0)
142 call readi(controlcard,'IPRINT',iprint,0)
143 call readi(controlcard,'MAXGEN',maxgen,10000)
144 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
145 call readi(controlcard,"KDIAG",kdiag,0)
146 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
147 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
148 & write (iout,*) "RESCALE_MODE",rescale_mode
149 split_ene=index(controlcard,'SPLIT_ENE').gt.0
150 if (index(controlcard,'REGULAR').gt.0.0D0) then
151 call reada(controlcard,'WEIDIS',weidis,0.1D0)
155 if (index(controlcard,'CHECKGRAD').gt.0) then
157 if (index(controlcard,'CART').gt.0) then
159 elseif (index(controlcard,'CARINT').gt.0) then
164 elseif (index(controlcard,'THREAD').gt.0) then
166 call readi(controlcard,'THREAD',nthread,0)
167 if (nthread.gt.0) then
168 call reada(controlcard,'WEIDIS',weidis,0.1D0)
171 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
172 stop 'Error termination in Read_Control.'
174 else if (index(controlcard,'MCMA').gt.0) then
176 else if (index(controlcard,'MCEE').gt.0) then
178 else if (index(controlcard,'MULTCONF').gt.0) then
180 else if (index(controlcard,'MAP').gt.0) then
182 call readi(controlcard,'MAP',nmap,0)
183 else if (index(controlcard,'CSA').gt.0) then
185 crc else if (index(controlcard,'ZSCORE').gt.0) then
187 crc ZSCORE is rm from UNRES, modecalc=9 is available
190 cfcm else if (index(controlcard,'MCMF').gt.0) then
192 else if (index(controlcard,'SOFTREG').gt.0) then
194 else if (index(controlcard,'CHECK_BOND').gt.0) then
196 else if (index(controlcard,'TEST').gt.0) then
198 else if (index(controlcard,'MD').gt.0) then
200 else if (index(controlcard,'RE ').gt.0) then
204 lmuca=index(controlcard,'MUCA').gt.0
205 call readi(controlcard,'MUCADYN',mucadyn,0)
206 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
207 if (lmuca .and. (me.eq.king .or. .not.out1file ))
209 write (iout,*) 'MUCADYN=',mucadyn
210 write (iout,*) 'MUCASMOOTH=',muca_smooth
213 iscode=index(controlcard,'ONE_LETTER')
214 indphi=index(controlcard,'PHI')
215 indback=index(controlcard,'BACK')
216 iranconf=index(controlcard,'RAND_CONF')
217 i2ndstr=index(controlcard,'USE_SEC_PRED')
218 gradout=index(controlcard,'GRADOUT').gt.0
219 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
220 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
221 if (me.eq.king .or. .not.out1file )
222 & write (iout,*) "DISTCHAINMAX",distchainmax
224 if(me.eq.king.or..not.out1file)
225 & write (iout,'(2a)') diagmeth(kdiag),
226 & ' routine used to diagonalize matrices.'
229 c--------------------------------------------------------------------------
230 subroutine read_REMDpar
234 implicit real*8 (a-h,o-z)
236 include 'COMMON.IOUNITS'
237 include 'COMMON.TIME1'
240 include 'COMMON.LANGEVIN'
242 include 'COMMON.LANGEVIN.lang0'
244 include 'COMMON.INTERACT'
245 include 'COMMON.NAMES'
247 include 'COMMON.REMD'
248 include 'COMMON.CONTROL'
249 include 'COMMON.SETUP'
251 character*320 controlcard
252 character*3200 controlcard1
253 integer iremd_m_total
255 if(me.eq.king.or..not.out1file)
256 & write (iout,*) "REMD setup"
258 call card_concat(controlcard)
259 call readi(controlcard,"NREP",nrep,3)
260 call readi(controlcard,"NSTEX",nstex,1000)
261 call reada(controlcard,"RETMIN",retmin,10.0d0)
262 call reada(controlcard,"RETMAX",retmax,1000.0d0)
263 mremdsync=(index(controlcard,'SYNC').gt.0)
264 call readi(controlcard,"NSYN",i_sync_step,100)
265 restart1file=(index(controlcard,'REST1FILE').gt.0)
266 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
267 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
268 if(max_cache_traj_use.gt.max_cache_traj)
269 & max_cache_traj_use=max_cache_traj
270 if(me.eq.king.or..not.out1file) then
271 cd if (traj1file) then
272 crc caching is in testing - NTWX is not ignored
273 cd write (iout,*) "NTWX value is ignored"
274 cd write (iout,*) " trajectory is stored to one file by master"
275 cd write (iout,*) " before exchange at NSTEX intervals"
277 write (iout,*) "NREP= ",nrep
278 write (iout,*) "NSTEX= ",nstex
279 write (iout,*) "SYNC= ",mremdsync
280 write (iout,*) "NSYN= ",i_sync_step
281 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
284 if (index(controlcard,'TLIST').gt.0) then
286 call card_concat(controlcard1)
287 read(controlcard1,*) (remd_t(i),i=1,nrep)
288 if(me.eq.king.or..not.out1file)
289 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
292 if (index(controlcard,'MLIST').gt.0) then
294 call card_concat(controlcard1)
295 read(controlcard1,*) (remd_m(i),i=1,nrep)
296 if(me.eq.king.or..not.out1file) then
297 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
300 iremd_m_total=iremd_m_total+remd_m(i)
302 write (iout,*) 'Total number of replicas ',iremd_m_total
305 if(me.eq.king.or..not.out1file)
306 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
309 c--------------------------------------------------------------------------
310 subroutine read_MDpar
314 implicit real*8 (a-h,o-z)
316 include 'COMMON.IOUNITS'
317 include 'COMMON.TIME1'
320 include 'COMMON.LANGEVIN'
322 include 'COMMON.LANGEVIN.lang0'
324 include 'COMMON.INTERACT'
325 include 'COMMON.NAMES'
327 include 'COMMON.SETUP'
328 include 'COMMON.CONTROL'
329 include 'COMMON.SPLITELE'
331 character*320 controlcard
333 call card_concat(controlcard)
334 call readi(controlcard,"NSTEP",n_timestep,1000000)
335 call readi(controlcard,"NTWE",ntwe,100)
336 call readi(controlcard,"NTWX",ntwx,1000)
337 call reada(controlcard,"DT",d_time,1.0d-1)
338 call reada(controlcard,"DVMAX",dvmax,2.0d1)
339 call reada(controlcard,"DAMAX",damax,1.0d1)
340 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
341 call readi(controlcard,"LANG",lang,0)
342 RESPA = index(controlcard,"RESPA") .gt. 0
343 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
344 ntime_split0=ntime_split
345 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
346 ntime_split0=ntime_split
347 call reada(controlcard,"R_CUT",r_cut,2.0d0)
348 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
349 rest = index(controlcard,"REST").gt.0
350 tbf = index(controlcard,"TBF").gt.0
351 usampl = index(controlcard,"USAMPL").gt.0
353 mdpdb = index(controlcard,"MDPDB").gt.0
354 call reada(controlcard,"T_BATH",t_bath,300.0d0)
355 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
356 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
357 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
358 if (count_reset_moment.eq.0) count_reset_moment=1000000000
359 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
360 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
361 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
362 if (count_reset_vel.eq.0) count_reset_vel=1000000000
363 large = index(controlcard,"LARGE").gt.0
364 print_compon = index(controlcard,"PRINT_COMPON").gt.0
365 rattle = index(controlcard,"RATTLE").gt.0
366 c if performing umbrella sampling, fragments constrained are read from the fragment file
372 if(me.eq.king.or..not.out1file) then
374 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
376 write (iout,'(a)') "The units are:"
377 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
378 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
379 & " acceleration: angstrom/(48.9 fs)**2"
380 write (iout,'(a)') "energy: kcal/mol, temperature: K"
382 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
383 write (iout,'(a60,f10.5,a)')
384 & "Initial time step of numerical integration:",d_time,
386 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
388 write (iout,'(2a,i4,a)')
389 & "A-MTS algorithm used; initial time step for fast-varying",
390 & " short-range forces split into",ntime_split," steps."
391 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
392 & r_cut," lambda",rlamb
394 write (iout,'(2a,f10.5)')
395 & "Maximum acceleration threshold to reduce the time step",
396 & "/increase split number:",damax
397 write (iout,'(2a,f10.5)')
398 & "Maximum predicted energy drift to reduce the timestep",
399 & "/increase split number:",edriftmax
400 write (iout,'(a60,f10.5)')
401 & "Maximum velocity threshold to reduce velocities:",dvmax
402 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
403 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
404 if (rattle) write (iout,'(a60)')
405 & "Rattle algorithm used to constrain the virtual bonds"
409 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
410 call reada(controlcard,"RWAT",rwat,1.4d0)
411 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
412 surfarea=index(controlcard,"SURFAREA").gt.0
413 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
414 if(me.eq.king.or..not.out1file)then
415 write (iout,'(/a,$)') "Langevin dynamics calculation"
418 & " with direct integration of Langevin equations"
419 else if (lang.eq.2) then
420 write (iout,'(a/)') " with TINKER stochasic MD integrator"
421 else if (lang.eq.3) then
422 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
423 else if (lang.eq.4) then
424 write (iout,'(a/)') " in overdamped mode"
426 write (iout,'(//a,i5)')
427 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
430 write (iout,'(a60,f10.5)') "Temperature:",t_bath
431 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
432 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
433 write (iout,'(a60,f10.5)')
434 & "Scaling factor of the friction forces:",scal_fric
435 if (surfarea) write (iout,'(2a,i10,a)')
436 & "Friction coefficients will be scaled by solvent-accessible",
437 & " surface area every",reset_fricmat," steps."
439 c Calculate friction coefficients and bounds of stochastic forces
440 eta=6*pi*cPoise*etawat
441 if(me.eq.king.or..not.out1file)
442 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
444 gamp=scal_fric*(pstok+rwat)*eta
445 stdfp=dsqrt(2*Rb*t_bath/d_time)
447 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
448 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
450 if(me.eq.king.or..not.out1file)then
451 write (iout,'(/2a/)')
452 & "Radii of site types and friction coefficients and std's of",
453 & " stochastic forces of fully exposed sites"
454 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
456 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
457 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
461 if(me.eq.king.or..not.out1file)then
462 write (iout,'(a)') "Berendsen bath calculation"
463 write (iout,'(a60,f10.5)') "Temperature:",t_bath
464 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
466 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
467 & count_reset_moment," steps"
469 & write (iout,'(a,i10,a)')
470 & "Velocities will be reset at random every",count_reset_vel,
474 if(me.eq.king.or..not.out1file)
475 & write (iout,'(a31)') "Microcanonical mode calculation"
477 if(me.eq.king.or..not.out1file)then
478 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
480 write(iout,*) "MD running with constraints."
481 write(iout,*) "Equilibration time ", eq_time, " mtus."
482 write(iout,*) "Constraining ", nfrag," fragments."
483 write(iout,*) "Length of each fragment, weight and q0:"
485 write (iout,*) "Set of restraints #",iset
487 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
488 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
490 write(iout,*) "constraints between ", npair, "fragments."
491 write(iout,*) "constraint pairs, weights and q0:"
493 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
494 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
496 write(iout,*) "angle constraints within ", nfrag_back,
497 & "backbone fragments."
498 write(iout,*) "fragment, weights:"
500 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
501 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
502 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
505 iset=mod(kolor,nset)+1
508 if(me.eq.king.or..not.out1file)
509 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
512 c------------------------------------------------------------------------------
515 C Read molecular data.
517 implicit real*8 (a-h,o-z)
523 include 'COMMON.IOUNITS'
526 include 'COMMON.INTERACT'
527 include 'COMMON.LOCAL'
528 include 'COMMON.NAMES'
529 include 'COMMON.CHAIN'
530 include 'COMMON.FFIELD'
531 include 'COMMON.SBRIDGE'
532 include 'COMMON.HEADER'
533 include 'COMMON.CONTROL'
534 include 'COMMON.DBASE'
535 include 'COMMON.THREAD'
536 include 'COMMON.CONTACTS'
537 include 'COMMON.TORCNSTR'
538 include 'COMMON.TIME1'
539 include 'COMMON.BOUNDS'
541 include 'COMMON.SETUP'
542 character*4 sequence(maxres)
544 double precision x(maxvar)
545 character*256 pdbfile
546 character*400 weightcard
547 character*80 weightcard_t,ucase
548 dimension itype_pdb(maxres)
549 common /pizda/ itype_pdb
550 logical seq_comp,fail
551 double precision energia(0:n_ene)
557 C Read weights of the subsequent energy terms.
558 call card_concat(weightcard)
559 call reada(weightcard,'WLONG',wlong,1.0D0)
560 call reada(weightcard,'WSC',wsc,wlong)
561 call reada(weightcard,'WSCP',wscp,wlong)
562 call reada(weightcard,'WELEC',welec,1.0D0)
563 call reada(weightcard,'WVDWPP',wvdwpp,welec)
564 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
565 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
566 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
567 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
568 call reada(weightcard,'WTURN3',wturn3,1.0D0)
569 call reada(weightcard,'WTURN4',wturn4,1.0D0)
570 call reada(weightcard,'WTURN6',wturn6,1.0D0)
571 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
572 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
573 call reada(weightcard,'WBOND',wbond,1.0D0)
574 call reada(weightcard,'WTOR',wtor,1.0D0)
575 call reada(weightcard,'WTORD',wtor_d,1.0D0)
576 call reada(weightcard,'WANG',wang,1.0D0)
577 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
578 call reada(weightcard,'SCAL14',scal14,0.4D0)
579 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
580 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
581 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
582 call reada(weightcard,'TEMP0',temp0,300.0d0)
583 if (index(weightcard,'SOFT').gt.0) ipot=6
584 C 12/1/95 Added weight for the multi-body term WCORR
585 call reada(weightcard,'WCORRH',wcorr,1.0D0)
586 if (wcorr4.gt.0.0d0) wcorr=wcorr4
606 if(me.eq.king.or..not.out1file)
607 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
608 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
610 10 format (/'Energy-term weights (unscaled):'//
611 & 'WSCC= ',f10.6,' (SC-SC)'/
612 & 'WSCP= ',f10.6,' (SC-p)'/
613 & 'WELEC= ',f10.6,' (p-p electr)'/
614 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
615 & 'WBOND= ',f10.6,' (stretching)'/
616 & 'WANG= ',f10.6,' (bending)'/
617 & 'WSCLOC= ',f10.6,' (SC local)'/
618 & 'WTOR= ',f10.6,' (torsional)'/
619 & 'WTORD= ',f10.6,' (double torsional)'/
620 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
621 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
622 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
623 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
624 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
625 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
626 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
627 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
628 & 'WTURN6= ',f10.6,' (turns, 6th order)')
629 if(me.eq.king.or..not.out1file)then
630 if (wcorr4.gt.0.0d0) then
631 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
632 & 'between contact pairs of peptide groups'
633 write (iout,'(2(a,f5.3/))')
634 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
635 & 'Range of quenching the correlation terms:',2*delt_corr
636 else if (wcorr.gt.0.0d0) then
637 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
638 & 'between contact pairs of peptide groups'
640 write (iout,'(a,f8.3)')
641 & 'Scaling factor of 1,4 SC-p interactions:',scal14
642 write (iout,'(a,f8.3)')
643 & 'General scaling factor of SC-p interactions:',scalscp
645 r0_corr=cutoff_corr-delt_corr
647 aad(i,1)=scalscp*aad(i,1)
648 aad(i,2)=scalscp*aad(i,2)
649 bad(i,1)=scalscp*bad(i,1)
650 bad(i,2)=scalscp*bad(i,2)
652 call rescale_weights(t_bath)
653 if(me.eq.king.or..not.out1file)
654 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
655 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
657 22 format (/'Energy-term weights (scaled):'//
658 & 'WSCC= ',f10.6,' (SC-SC)'/
659 & 'WSCP= ',f10.6,' (SC-p)'/
660 & 'WELEC= ',f10.6,' (p-p electr)'/
661 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
662 & 'WBOND= ',f10.6,' (stretching)'/
663 & 'WANG= ',f10.6,' (bending)'/
664 & 'WSCLOC= ',f10.6,' (SC local)'/
665 & 'WTOR= ',f10.6,' (torsional)'/
666 & 'WTORD= ',f10.6,' (double torsional)'/
667 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
668 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
669 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
670 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
671 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
672 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
673 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
674 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
675 & 'WTURN6= ',f10.6,' (turns, 6th order)')
676 if(me.eq.king.or..not.out1file)
677 & write (iout,*) "Reference temperature for weights calculation:",
679 call reada(weightcard,"D0CM",d0cm,3.78d0)
680 call reada(weightcard,"AKCM",akcm,15.1d0)
681 call reada(weightcard,"AKTH",akth,11.0d0)
682 call reada(weightcard,"AKCT",akct,12.0d0)
683 call reada(weightcard,"V1SS",v1ss,-1.08d0)
684 call reada(weightcard,"V2SS",v2ss,7.61d0)
685 call reada(weightcard,"V3SS",v3ss,13.7d0)
686 call reada(weightcard,"EBR",ebr,-5.50D0)
687 call reada(weightcard,"ATRISS",atriss,0.301D0)
688 call reada(weightcard,"BTRISS",btriss,0.021D0)
689 call reada(weightcard,"CTRISS",ctriss,1.001D0)
690 call reada(weightcard,"DTRISS",dtriss,1.001D0)
691 write (iout,*) "ATRISS=", atriss
692 write (iout,*) "BTRISS=", btriss
693 write (iout,*) "CTRISS=", ctriss
694 write (iout,*) "DTRISS=", dtriss
695 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
697 dyn_ss_mask(i)=.false.
701 dyn_ssbond_ij(i,j)=1.0d300
704 call reada(weightcard,"HT",Ht,0.0D0)
706 ss_depth=ebr/wsc-0.25*eps(1,1)
707 Ht=Ht/wsc-0.25*eps(1,1)
708 akcm=akcm*wstrain/wsc
709 akth=akth*wstrain/wsc
710 akct=akct*wstrain/wsc
711 v1ss=v1ss*wstrain/wsc
712 v2ss=v2ss*wstrain/wsc
713 v3ss=v3ss*wstrain/wsc
715 if (wstrain.ne.0.0) then
716 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
722 if(me.eq.king.or..not.out1file) then
723 write (iout,*) "Parameters of the SS-bond potential:"
724 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
726 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
727 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
728 write (iout,*)" HT",Ht
729 print *,'indpdb=',indpdb,' pdbref=',pdbref
731 if (indpdb.gt.0 .or. pdbref) then
732 read(inp,'(a)') pdbfile
733 if(me.eq.king.or..not.out1file)
734 & write (iout,'(2a)') 'PDB data will be read from file ',
735 & pdbfile(:ilen(pdbfile))
736 open(ipdbin,file=pdbfile,status='old',err=33)
738 33 write (iout,'(a)') 'Error opening PDB file.'
741 c write (iout,*) 'Begin reading pdb data'
744 c write (iout,*) 'Finished reading pdb data'
746 if(me.eq.king.or..not.out1file)
747 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
748 & ' nstart_sup=',nstart_sup
750 itype_pdb(i)=itype(i)
754 nct=nstart_sup+nsup-1
755 call contact(.false.,ncont_ref,icont_ref,co)
758 if(me.eq.king.or..not.out1file)
759 & write(iout,*)'Adding sidechains'
763 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
766 do while (fail.and.nsi.le.maxsi)
767 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
770 if(fail) write(iout,*)'Adding sidechain failed for res ',
771 & i,' after ',nsi,' trials'
776 if (indpdb.eq.0) then
777 C Read sequence if not taken from the pdb file.
779 c print *,'nres=',nres
780 if (iscode.gt.0) then
781 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
783 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
785 C Convert sequence to numeric code
787 itype(i)=rescode(i,sequence(i),iscode)
789 C Assign initial virtual bond lengths
795 vbld(i+nres)=dsc(iabs(itype(i)))
796 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
797 c write (iout,*) "i",i," itype",itype(i),
798 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
802 c print '(20i4)',(itype(i),i=1,nres)
805 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
807 if (itype(i).eq.ntyp1) then
811 else if (iabs(itype(i+1)).ne.20) then
813 else if (iabs(itype(i)).ne.20) then
820 if(me.eq.king.or..not.out1file)then
821 write (iout,*) "ITEL"
823 write (iout,*) i,itype(i),itel(i)
825 print *,'Call Read_Bridge.'
828 C 8/13/98 Set limits to generating the dihedral angles
833 read (inp,*) ndih_constr
834 if (ndih_constr.gt.0) then
836 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
838 if(me.eq.king.or..not.out1file)then
840 & 'There are',ndih_constr,' constraints on phi angles.'
842 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
847 phi0(i)=deg2rad*phi0(i)
848 drange(i)=deg2rad*drange(i)
850 C if(me.eq.king.or..not.out1file)
851 C & write (iout,*) 'FTORS',ftors
854 phibound(1,ii) = phi0(i)-drange(i)
855 phibound(2,ii) = phi0(i)+drange(i)
858 C first setting the theta boundaries to 0 to pi
859 C this mean that there is no energy penalty for any angle occuring
864 C begin reading theta constrains this is quartic constrains allowing to
865 C have smooth second derivative
866 if (with_theta_constr) then
867 C with_theta_constr is keyword allowing for occurance of theta constrains
868 read (inp,*) ntheta_constr
869 C ntheta_constr is the number of theta constrains
870 if (ntheta_constr.gt.0) then
872 read (inp,*) (itheta_constr(i),theta_constr0(i),
873 & theta_drange(i),for_thet_constr(i),
875 C the above code reads from 1 to ntheta_constr
876 C itheta_constr(i) residue i for which is theta_constr
877 C theta_constr0 the global minimum value
878 C theta_drange is range for which there is no energy penalty
879 C for_thet_constr is the force constant for quartic energy penalty
881 if(me.eq.king.or..not.out1file)then
883 & 'There are',ntheta_constr,' constraints on phi angles.'
885 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
891 theta_constr0(i)=deg2rad*theta_constr0(i)
892 theta_drange(i)=deg2rad*theta_drange(i)
894 C if(me.eq.king.or..not.out1file)
895 C & write (iout,*) 'FTORS',ftors
897 ii = itheta_constr(i)
898 thetabound(1,ii) = phi0(i)-drange(i)
899 thetabound(2,ii) = phi0(i)+drange(i)
901 endif ! ntheta_constr.gt.0
902 endif! with_theta_constr
904 C with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
905 C write (iout,*) "with_dihed_constr ",with_dihed_constr
910 write (iout,'(a)') 'Boundaries in phi angle sampling:'
912 write (iout,'(a3,i5,2f10.1)')
913 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
919 cd print *,'NNT=',NNT,' NCT=',NCT
920 if (itype(1).eq.ntyp1) nnt=2
921 if (itype(nres).eq.ntyp1) nct=nct-1
923 if(me.eq.king.or..not.out1file)
924 & write (iout,'(a,i3)') 'nsup=',nsup
926 if (nsup.le.(nct-nnt+1)) then
927 do i=0,nct-nnt+1-nsup
928 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
934 & 'Error - sequences to be superposed do not match.'
937 do i=0,nsup-(nct-nnt+1)
938 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
940 nstart_sup=nstart_sup+i
946 & 'Error - sequences to be superposed do not match.'
949 if (nsup.eq.0) nsup=nct-nnt
950 if (nstart_sup.eq.0) nstart_sup=nnt
951 if (nstart_seq.eq.0) nstart_seq=nnt
952 if(me.eq.king.or..not.out1file)
953 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
954 & ' nstart_seq=',nstart_seq
956 c--- Zscore rms -------
957 if (nz_start.eq.0) nz_start=nnt
958 if (nz_end.eq.0 .and. nsup.gt.0) then
960 else if (nz_end.eq.0) then
963 if(me.eq.king.or..not.out1file)then
964 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
965 write (iout,*) 'IZ_SC=',iz_sc
967 c----------------------
970 if (.not.pdbref) then
971 call read_angles(inp,*38)
973 38 write (iout,'(a)') 'Error reading reference structure.'
975 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
976 stop 'Error reading reference structure'
980 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
990 call contact(.true.,ncont_ref,icont_ref,co)
994 write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
996 if (constr_dist.gt.0) call read_dist_constr
997 write (iout,*) "After read_dist_constr nhpb",nhpb
999 if(me.eq.king.or..not.out1file)
1000 & write (iout,*) 'Contact order:',co
1002 if(me.eq.king.or..not.out1file)
1003 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1006 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1008 if(me.eq.king.or..not.out1file)
1009 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1010 & icont_ref(1,i),' ',
1011 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1015 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1016 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1017 & modecalc.ne.10) then
1018 C If input structure hasn't been supplied from the PDB file read or generate
1020 if (iranconf.eq.0 .and. .not. extconf) then
1021 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1022 & write (iout,'(a)') 'Initial geometry will be read in.'
1024 read(inp,'(8f10.5)',end=36,err=36)
1025 & ((c(l,k),l=1,3),k=1,nres),
1026 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1027 write (iout,*) "Exit READ_CART"
1028 write (iout,'(8f10.5)')
1029 & ((c(l,k),l=1,3),k=1,nres),
1030 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1031 call int_from_cart1(.true.)
1032 write (iout,*) "Finish INT_TO_CART"
1035 dc(j,i)=c(j,i+1)-c(j,i)
1036 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1040 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1042 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1043 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1049 call read_angles(inp,*36)
1052 36 write (iout,'(a)') 'Error reading angle file.'
1054 call mpi_finalize( MPI_COMM_WORLD,IERR )
1056 stop 'Error reading angle file.'
1058 else if (extconf) then
1059 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1060 & write (iout,'(a)') 'Extended chain initial geometry.'
1062 theta(i)=90d0*deg2rad
1065 phi(i)=180d0*deg2rad
1068 alph(i)=110d0*deg2rad
1071 omeg(i)=-120d0*deg2rad
1072 if (itype(i).le.0) omeg(i)=-omeg(i)
1075 if(me.eq.king.or..not.out1file)
1076 & write (iout,'(a)') 'Random-generated initial geometry.'
1080 if (me.eq.king .or. fg_rank.eq.0 .and. (
1081 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1085 call gen_rand_conf(itmp,*30)
1087 30 write (iout,*) 'Failed to generate random conformation',
1088 & ', itrial=',itrial
1089 write (*,*) 'Processor:',me,
1090 & ' Failed to generate random conformation',
1100 write (iout,'(a,i3,a)') 'Processor:',me,
1101 & ' error in generating random conformation.'
1102 write (*,'(a,i3,a)') 'Processor:',me,
1103 & ' error in generating random conformation.'
1106 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1112 call gen_rand_conf(itmp,*30)
1114 30 write (iout,*) 'Failed to generate random conformation',
1115 & ', itrial=',itrial
1116 write (*,*) 'Failed to generate random conformation',
1117 & ', itrial=',itrial
1119 write (iout,'(a,i3,a)') 'Processor:',me,
1120 & ' error in generating random conformation.'
1121 write (*,'(a,i3,a)') 'Processor:',me,
1122 & ' error in generating random conformation.'
1127 elseif (modecalc.eq.4) then
1128 read (inp,'(a)') intinname
1129 open (intin,file=intinname,status='old',err=333)
1130 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1131 & write (iout,'(a)') 'intinname',intinname
1132 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1134 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1136 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1138 stop 'Error opening angle file.'
1142 C Generate distance constraints, if the PDB structure is to be regularized.
1143 if (nthread.gt.0) then
1144 call read_threadbase
1147 if (me.eq.king .or. .not. out1file)
1149 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1150 write (iout,'(/a,i3,a)')
1151 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1152 write (iout,'(20i4)') (iss(i),i=1,ns)
1154 write(iout,*)"Running with dynamic disulfide-bond formation"
1156 write (iout,'(/a/)') 'Pre-formed links are:'
1162 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1163 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1169 if (ns.gt.0.and.dyn_ss) then
1173 forcon(i-nss)=forcon(i)
1180 dyn_ss_mask(iss(i))=.true.
1183 if (i2ndstr.gt.0) call secstrp2dihc
1184 c call geom_to_var(nvar,x)
1185 c call etotal(energia(0))
1186 c call enerprint(energia(0))
1187 c call briefout(0,etot)
1189 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1190 cd write (iout,'(a)') 'Variable list:'
1191 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1193 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1194 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1195 & 'Processor',myrank,': end reading molecular data.'
1200 c--------------------------------------------------------------------------
1201 logical function seq_comp(itypea,itypeb,length)
1203 integer length,itypea(length),itypeb(length)
1206 if (itypea(i).ne.itypeb(i)) then
1214 c-----------------------------------------------------------------------------
1215 subroutine read_bridge
1216 C Read information about disulfide bridges.
1217 implicit real*8 (a-h,o-z)
1218 include 'DIMENSIONS'
1222 include 'COMMON.IOUNITS'
1223 include 'COMMON.GEO'
1224 include 'COMMON.VAR'
1225 include 'COMMON.INTERACT'
1226 include 'COMMON.LOCAL'
1227 include 'COMMON.NAMES'
1228 include 'COMMON.CHAIN'
1229 include 'COMMON.FFIELD'
1230 include 'COMMON.SBRIDGE'
1231 include 'COMMON.HEADER'
1232 include 'COMMON.CONTROL'
1233 include 'COMMON.DBASE'
1234 include 'COMMON.THREAD'
1235 include 'COMMON.TIME1'
1236 include 'COMMON.SETUP'
1237 C Read bridging residues.
1238 read (inp,*) ns,(iss(i),i=1,ns)
1240 if(me.eq.king.or..not.out1file)
1241 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1242 C Check whether the specified bridging residues are cystines.
1244 if (itype(iss(i)).ne.1) then
1245 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1246 & 'Do you REALLY think that the residue ',
1247 & restyp(itype(iss(i))),i,
1248 & ' can form a disulfide bridge?!!!'
1249 write (*,'(2a,i3,a)')
1250 & 'Do you REALLY think that the residue ',
1251 & restyp(itype(iss(i))),i,
1252 & ' can form a disulfide bridge?!!!'
1254 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1259 C Read preformed bridges.
1261 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1263 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1266 C Check if the residues involved in bridges are in the specified list of
1267 C bridging residues.
1270 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1271 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1272 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1273 & ' contains residues present in other pairs.'
1274 write (*,'(a,i3,a)') 'Disulfide pair',i,
1275 & ' contains residues present in other pairs.'
1277 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1283 if (ihpb(i).eq.iss(j)) goto 10
1285 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1288 if (jhpb(i).eq.iss(j)) goto 20
1290 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1296 ihpb(i)=ihpb(i)+nres
1297 jhpb(i)=jhpb(i)+nres
1303 c----------------------------------------------------------------------------
1304 subroutine read_x(kanal,*)
1305 implicit real*8 (a-h,o-z)
1306 include 'DIMENSIONS'
1307 include 'COMMON.GEO'
1308 include 'COMMON.VAR'
1309 include 'COMMON.CHAIN'
1310 include 'COMMON.IOUNITS'
1311 include 'COMMON.CONTROL'
1312 include 'COMMON.LOCAL'
1313 include 'COMMON.INTERACT'
1314 c Read coordinates from input
1316 read(kanal,'(8f10.5)',end=10,err=10)
1317 & ((c(l,k),l=1,3),k=1,nres),
1318 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1321 c(j,2*nres)=c(j,nres)
1323 call int_from_cart1(.false.)
1326 dc(j,i)=c(j,i+1)-c(j,i)
1327 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1331 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1333 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1334 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1342 c----------------------------------------------------------------------------
1343 subroutine read_threadbase
1344 implicit real*8 (a-h,o-z)
1345 include 'DIMENSIONS'
1346 include 'COMMON.IOUNITS'
1347 include 'COMMON.GEO'
1348 include 'COMMON.VAR'
1349 include 'COMMON.INTERACT'
1350 include 'COMMON.LOCAL'
1351 include 'COMMON.NAMES'
1352 include 'COMMON.CHAIN'
1353 include 'COMMON.FFIELD'
1354 include 'COMMON.SBRIDGE'
1355 include 'COMMON.HEADER'
1356 include 'COMMON.CONTROL'
1357 include 'COMMON.DBASE'
1358 include 'COMMON.THREAD'
1359 include 'COMMON.TIME1'
1360 C Read pattern database for threading.
1361 read (icbase,*) nseq
1363 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1364 & nres_base(2,i),nres_base(3,i)
1365 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1367 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1368 c & nres_base(2,i),nres_base(3,i)
1369 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1373 if (weidis.eq.0.0D0) weidis=0.1D0
1382 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1383 write (iout,'(a,i5)') 'nexcl: ',nexcl
1384 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1387 c------------------------------------------------------------------------------
1388 subroutine setup_var
1389 implicit real*8 (a-h,o-z)
1390 include 'DIMENSIONS'
1391 include 'COMMON.IOUNITS'
1392 include 'COMMON.GEO'
1393 include 'COMMON.VAR'
1394 include 'COMMON.INTERACT'
1395 include 'COMMON.LOCAL'
1396 include 'COMMON.NAMES'
1397 include 'COMMON.CHAIN'
1398 include 'COMMON.FFIELD'
1399 include 'COMMON.SBRIDGE'
1400 include 'COMMON.HEADER'
1401 include 'COMMON.CONTROL'
1402 include 'COMMON.DBASE'
1403 include 'COMMON.THREAD'
1404 include 'COMMON.TIME1'
1405 C Set up variable list.
1411 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1413 ialph(i,1)=nvar+nside
1417 if (indphi.gt.0) then
1419 else if (indback.gt.0) then
1424 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1427 c----------------------------------------------------------------------------
1428 subroutine gen_dist_constr
1429 C Generate CA distance constraints.
1430 implicit real*8 (a-h,o-z)
1431 include 'DIMENSIONS'
1432 include 'COMMON.IOUNITS'
1433 include 'COMMON.GEO'
1434 include 'COMMON.VAR'
1435 include 'COMMON.INTERACT'
1436 include 'COMMON.LOCAL'
1437 include 'COMMON.NAMES'
1438 include 'COMMON.CHAIN'
1439 include 'COMMON.FFIELD'
1440 include 'COMMON.SBRIDGE'
1441 include 'COMMON.HEADER'
1442 include 'COMMON.CONTROL'
1443 include 'COMMON.DBASE'
1444 include 'COMMON.THREAD'
1445 include 'COMMON.TIME1'
1446 dimension itype_pdb(maxres)
1447 common /pizda/ itype_pdb
1449 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1450 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1451 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1453 do i=nstart_sup,nstart_sup+nsup-1
1454 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1455 cd & ' seq_pdb', restyp(itype_pdb(i))
1456 do j=i+2,nstart_sup+nsup-1
1458 ihpb(nhpb)=i+nstart_seq-nstart_sup
1459 jhpb(nhpb)=j+nstart_seq-nstart_sup
1461 dhpb(nhpb)=dist(i,j)
1464 cd write (iout,'(a)') 'Distance constraints:'
1469 cd if (ii.gt.nres) then
1474 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1475 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1476 cd & dhpb(i),forcon(i)
1480 c----------------------------------------------------------------------------
1482 implicit real*8 (a-h,o-z)
1483 include 'DIMENSIONS'
1484 include 'COMMON.MAP'
1485 include 'COMMON.IOUNITS'
1486 character*3 angid(4) /'THE','PHI','ALP','OME'/
1487 character*80 mapcard,ucase
1489 read (inp,'(a)') mapcard
1490 mapcard=ucase(mapcard)
1491 if (index(mapcard,'PHI').gt.0) then
1493 else if (index(mapcard,'THE').gt.0) then
1495 else if (index(mapcard,'ALP').gt.0) then
1497 else if (index(mapcard,'OME').gt.0) then
1500 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1501 stop 'Error - illegal variable spec in MAP card.'
1503 call readi (mapcard,'RES1',res1(imap),0)
1504 call readi (mapcard,'RES2',res2(imap),0)
1505 if (res1(imap).eq.0) then
1506 res1(imap)=res2(imap)
1507 else if (res2(imap).eq.0) then
1508 res2(imap)=res1(imap)
1510 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1512 & 'Error - illegal definition of variable group in MAP.'
1513 stop 'Error - illegal definition of variable group in MAP.'
1515 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1516 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1517 call readi(mapcard,'NSTEP',nstep(imap),0)
1518 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1520 & 'Illegal boundary and/or step size specification in MAP.'
1521 stop 'Illegal boundary and/or step size specification in MAP.'
1526 c----------------------------------------------------------------------------
1528 implicit real*8 (a-h,o-z)
1529 include 'DIMENSIONS'
1530 include 'COMMON.IOUNITS'
1531 include 'COMMON.GEO'
1532 include 'COMMON.CSA'
1533 include 'COMMON.BANK'
1534 include 'COMMON.CONTROL'
1536 character*620 mcmcard
1537 call card_concat(mcmcard)
1539 call readi(mcmcard,'NCONF',nconf,50)
1540 call readi(mcmcard,'NADD',nadd,0)
1541 call readi(mcmcard,'JSTART',jstart,1)
1542 call readi(mcmcard,'JEND',jend,1)
1543 call readi(mcmcard,'NSTMAX',nstmax,500000)
1544 call readi(mcmcard,'N0',n0,1)
1545 call readi(mcmcard,'N1',n1,6)
1546 call readi(mcmcard,'N2',n2,4)
1547 call readi(mcmcard,'N3',n3,0)
1548 call readi(mcmcard,'N4',n4,0)
1549 call readi(mcmcard,'N5',n5,0)
1550 call readi(mcmcard,'N6',n6,10)
1551 call readi(mcmcard,'N7',n7,0)
1552 call readi(mcmcard,'N8',n8,0)
1553 call readi(mcmcard,'N9',n9,0)
1554 call readi(mcmcard,'N14',n14,0)
1555 call readi(mcmcard,'N15',n15,0)
1556 call readi(mcmcard,'N16',n16,0)
1557 call readi(mcmcard,'N17',n17,0)
1558 call readi(mcmcard,'N18',n18,0)
1560 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1562 call readi(mcmcard,'NDIFF',ndiff,2)
1563 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1564 call readi(mcmcard,'IS1',is1,1)
1565 call readi(mcmcard,'IS2',is2,8)
1566 call readi(mcmcard,'NRAN0',nran0,4)
1567 call readi(mcmcard,'NRAN1',nran1,2)
1568 call readi(mcmcard,'IRR',irr,1)
1569 call readi(mcmcard,'NSEED',nseed,20)
1570 call readi(mcmcard,'NTOTAL',ntotal,10000)
1571 call reada(mcmcard,'CUT1',cut1,2.0d0)
1572 call reada(mcmcard,'CUT2',cut2,5.0d0)
1573 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1574 call readi(mcmcard,'ICMAX',icmax,3)
1575 call readi(mcmcard,'IRESTART',irestart,0)
1576 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1579 call reada(mcmcard,'DELE',dele,20.0d0)
1580 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1581 call readi(mcmcard,'IREF',iref,0)
1582 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1583 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1584 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1585 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1586 write (iout,*) "NCONF_IN",nconf_in
1589 c----------------------------------------------------------------------------
1590 cfmc subroutine mcmfread
1591 cfmc implicit real*8 (a-h,o-z)
1592 cfmc include 'DIMENSIONS'
1593 cfmc include 'COMMON.MCMF'
1594 cfmc include 'COMMON.IOUNITS'
1595 cfmc include 'COMMON.GEO'
1596 cfmc character*80 ucase
1597 cfmc character*620 mcmcard
1598 cfmc call card_concat(mcmcard)
1600 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1601 cfmc write(iout,*)'MAXRANT=',maxrant
1602 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1603 cfmc write(iout,*)'MAXFAM=',maxfam
1604 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1605 cfmc write(iout,*)'NNET1=',nnet1
1606 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1607 cfmc write(iout,*)'NNET2=',nnet2
1608 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1609 cfmc write(iout,*)'NNET3=',nnet3
1610 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1611 cfmc write(iout,*)'ILASTT=',ilastt
1612 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1613 cfmc write(iout,*)'MAXSTR=',maxstr
1614 cfmc maxstr_f=maxstr/maxfam
1615 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1616 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1617 cfmc write(iout,*)'NMCMF=',nmcmf
1618 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1619 cfmc write(iout,*)'IFOCUS=',ifocus
1620 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1621 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1622 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1623 cfmc write(iout,*)'INTPRT=',intprt
1624 cfmc call readi(mcmcard,'IPRT',iprt,100)
1625 cfmc write(iout,*)'IPRT=',iprt
1626 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1627 cfmc write(iout,*)'IMAXTR=',imaxtr
1628 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1629 cfmc write(iout,*)'MAXEVEN=',maxeven
1630 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1631 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1632 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1633 cfmc write(iout,*)'INIMIN=',inimin
1634 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1635 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1636 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1637 cfmc write(iout,*)'NTHREAD=',nthread
1638 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1639 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1640 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1641 cfmc write(iout,*)'MAXPERT=',maxpert
1642 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1643 cfmc write(iout,*)'IRMSD=',irmsd
1644 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1645 cfmc write(iout,*)'DENEMIN=',denemin
1646 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1647 cfmc write(iout,*)'RCUT1S=',rcut1s
1648 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1649 cfmc write(iout,*)'RCUT1E=',rcut1e
1650 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1651 cfmc write(iout,*)'RCUT2S=',rcut2s
1652 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1653 cfmc write(iout,*)'RCUT2E=',rcut2e
1654 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1655 cfmc write(iout,*)'DPERT1=',d_pert1
1656 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1657 cfmc write(iout,*)'DPERT1A=',d_pert1a
1658 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1659 cfmc write(iout,*)'DPERT2=',d_pert2
1660 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1661 cfmc write(iout,*)'DPERT2A=',d_pert2a
1662 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1663 cfmc write(iout,*)'DPERT2B=',d_pert2b
1664 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1665 cfmc write(iout,*)'DPERT2C=',d_pert2c
1666 cfmc d_pert1=deg2rad*d_pert1
1667 cfmc d_pert1a=deg2rad*d_pert1a
1668 cfmc d_pert2=deg2rad*d_pert2
1669 cfmc d_pert2a=deg2rad*d_pert2a
1670 cfmc d_pert2b=deg2rad*d_pert2b
1671 cfmc d_pert2c=deg2rad*d_pert2c
1672 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1673 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1674 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1675 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1676 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1677 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1678 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1679 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1680 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1681 cfmc write(iout,*)'RCUTINI=',rcutini
1682 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1683 cfmc write(iout,*)'GRAT=',grat
1684 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1685 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1689 c----------------------------------------------------------------------------
1691 implicit real*8 (a-h,o-z)
1692 include 'DIMENSIONS'
1693 include 'COMMON.MCM'
1694 include 'COMMON.MCE'
1695 include 'COMMON.IOUNITS'
1697 character*320 mcmcard
1698 call card_concat(mcmcard)
1699 call readi(mcmcard,'MAXACC',maxacc,100)
1700 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1701 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1702 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1703 call readi(mcmcard,'MAXREPM',maxrepm,200)
1704 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1705 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1706 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1707 call reada(mcmcard,'E_UP',e_up,5.0D0)
1708 call reada(mcmcard,'DELTE',delte,0.1D0)
1709 call readi(mcmcard,'NSWEEP',nsweep,5)
1710 call readi(mcmcard,'NSTEPH',nsteph,0)
1711 call readi(mcmcard,'NSTEPC',nstepc,0)
1712 call reada(mcmcard,'TMIN',tmin,298.0D0)
1713 call reada(mcmcard,'TMAX',tmax,298.0D0)
1714 call readi(mcmcard,'NWINDOW',nwindow,0)
1715 call readi(mcmcard,'PRINT_MC',print_mc,0)
1716 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1717 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1718 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1719 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1720 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1721 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1722 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1723 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1724 if (nwindow.gt.0) then
1725 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1727 winlen(i)=winend(i)-winstart(i)+1
1730 if (tmax.lt.tmin) tmax=tmin
1731 if (tmax.eq.tmin) then
1735 if (nstepc.gt.0 .and. nsteph.gt.0) then
1736 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1737 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1739 C Probabilities of different move types
1740 sumpro_type(0)=0.0D0
1741 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1742 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1743 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1744 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1745 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1746 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1747 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1749 print *,'i',i,' sumprotype',sumpro_type(i)
1750 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1751 print *,'i',i,' sumprotype',sumpro_type(i)
1755 c----------------------------------------------------------------------------
1756 subroutine read_minim
1757 implicit real*8 (a-h,o-z)
1758 include 'DIMENSIONS'
1759 include 'COMMON.MINIM'
1760 include 'COMMON.IOUNITS'
1762 character*320 minimcard
1763 call card_concat(minimcard)
1764 call readi(minimcard,'MAXMIN',maxmin,2000)
1765 call readi(minimcard,'MAXFUN',maxfun,5000)
1766 call readi(minimcard,'MINMIN',minmin,maxmin)
1767 call readi(minimcard,'MINFUN',minfun,maxmin)
1768 call reada(minimcard,'TOLF',tolf,1.0D-2)
1769 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1770 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1771 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1772 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1773 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1774 & 'Options in energy minimization:'
1775 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1776 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1777 & 'MinMin:',MinMin,' MinFun:',MinFun,
1778 & ' TolF:',TolF,' RTolF:',RTolF
1781 c----------------------------------------------------------------------------
1782 subroutine read_angles(kanal,*)
1783 implicit real*8 (a-h,o-z)
1784 include 'DIMENSIONS'
1785 include 'COMMON.GEO'
1786 include 'COMMON.VAR'
1787 include 'COMMON.CHAIN'
1788 include 'COMMON.IOUNITS'
1789 include 'COMMON.CONTROL'
1790 c Read angles from input
1792 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1793 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1794 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1795 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1798 c 9/7/01 avoid 180 deg valence angle
1799 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1801 theta(i)=deg2rad*theta(i)
1802 phi(i)=deg2rad*phi(i)
1803 alph(i)=deg2rad*alph(i)
1804 omeg(i)=deg2rad*omeg(i)
1809 c----------------------------------------------------------------------------
1810 subroutine reada(rekord,lancuch,wartosc,default)
1812 character*(*) rekord,lancuch
1813 double precision wartosc,default
1816 iread=index(rekord,lancuch)
1817 if (iread.eq.0) then
1821 iread=iread+ilen(lancuch)+1
1822 read (rekord(iread:),*,err=10,end=10) wartosc
1827 c----------------------------------------------------------------------------
1828 subroutine readi(rekord,lancuch,wartosc,default)
1830 character*(*) rekord,lancuch
1831 integer wartosc,default
1834 iread=index(rekord,lancuch)
1835 if (iread.eq.0) then
1839 iread=iread+ilen(lancuch)+1
1840 read (rekord(iread:),*,err=10,end=10) wartosc
1845 c----------------------------------------------------------------------------
1846 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1849 integer tablica(dim),default
1850 character*(*) rekord,lancuch
1857 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1858 if (iread.eq.0) return
1859 iread=iread+ilen(lancuch)+1
1860 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1863 c----------------------------------------------------------------------------
1864 subroutine multreada(rekord,lancuch,tablica,dim,default)
1867 double precision tablica(dim),default
1868 character*(*) rekord,lancuch
1875 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1876 if (iread.eq.0) return
1877 iread=iread+ilen(lancuch)+1
1878 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1881 c----------------------------------------------------------------------------
1882 subroutine openunits
1883 implicit real*8 (a-h,o-z)
1884 include 'DIMENSIONS'
1887 character*16 form,nodename
1890 include 'COMMON.SETUP'
1891 include 'COMMON.IOUNITS'
1893 include 'COMMON.CONTROL'
1894 integer lenpre,lenpot,ilen,lentmp
1896 character*3 out1file_text,ucase
1899 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1900 call getenv_loc("PREFIX",prefix)
1902 call getenv_loc("POT",pot)
1903 call getenv_loc("DIRTMP",tmpdir)
1904 call getenv_loc("CURDIR",curdir)
1905 call getenv_loc("OUT1FILE",out1file_text)
1906 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1907 out1file_text=ucase(out1file_text)
1908 if (out1file_text(1:1).eq."Y") then
1911 out1file=fg_rank.gt.0
1916 if (lentmp.gt.0) then
1917 write (*,'(80(1h!))')
1918 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1919 write (*,'(80(1h!))')
1920 write (*,*)"All output files will be on node /tmp directory."
1922 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1923 if (me.eq.king) then
1924 write (*,*) "The master node is ",nodename
1925 else if (fg_rank.eq.0) then
1926 write (*,*) "I am the CG slave node ",nodename
1928 write (*,*) "I am the FG slave node ",nodename
1931 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1932 lenpre = lentmp+lenpre+1
1934 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1935 C Get the names and open the input files
1936 #if defined(WINIFL) || defined(WINPGI)
1937 open(1,file=pref_orig(:ilen(pref_orig))//
1938 & '.inp',status='old',readonly,shared)
1939 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1940 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1941 C Get parameter filenames and open the parameter files.
1942 call getenv_loc('BONDPAR',bondname)
1943 open (ibond,file=bondname,status='old',readonly,shared)
1944 call getenv_loc('THETPAR',thetname)
1945 open (ithep,file=thetname,status='old',readonly,shared)
1946 call getenv_loc('ROTPAR',rotname)
1947 open (irotam,file=rotname,status='old',readonly,shared)
1948 call getenv_loc('TORPAR',torname)
1949 open (itorp,file=torname,status='old',readonly,shared)
1950 call getenv_loc('TORDPAR',tordname)
1951 open (itordp,file=tordname,status='old',readonly,shared)
1952 call getenv_loc('FOURIER',fouriername)
1953 open (ifourier,file=fouriername,status='old',readonly,shared)
1954 call getenv_loc('ELEPAR',elename)
1955 open (ielep,file=elename,status='old',readonly,shared)
1956 call getenv_loc('SIDEPAR',sidename)
1957 open (isidep,file=sidename,status='old',readonly,shared)
1958 #elif (defined CRAY) || (defined AIX)
1959 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1961 c print *,"Processor",myrank," opened file 1"
1962 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1963 c print *,"Processor",myrank," opened file 9"
1964 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1965 C Get parameter filenames and open the parameter files.
1966 call getenv_loc('BONDPAR',bondname)
1967 open (ibond,file=bondname,status='old',action='read')
1968 c print *,"Processor",myrank," opened file IBOND"
1969 call getenv_loc('THETPAR',thetname)
1970 open (ithep,file=thetname,status='old',action='read')
1971 c print *,"Processor",myrank," opened file ITHEP"
1972 call getenv_loc('ROTPAR',rotname)
1973 open (irotam,file=rotname,status='old',action='read')
1974 c print *,"Processor",myrank," opened file IROTAM"
1975 call getenv_loc('TORPAR',torname)
1976 open (itorp,file=torname,status='old',action='read')
1977 c print *,"Processor",myrank," opened file ITORP"
1978 call getenv_loc('TORDPAR',tordname)
1979 open (itordp,file=tordname,status='old',action='read')
1980 c print *,"Processor",myrank," opened file ITORDP"
1981 call getenv_loc('SCCORPAR',sccorname)
1982 open (isccor,file=sccorname,status='old',action='read')
1983 c print *,"Processor",myrank," opened file ISCCOR"
1984 call getenv_loc('FOURIER',fouriername)
1985 open (ifourier,file=fouriername,status='old',action='read')
1986 c print *,"Processor",myrank," opened file IFOURIER"
1987 call getenv_loc('ELEPAR',elename)
1988 open (ielep,file=elename,status='old',action='read')
1989 c print *,"Processor",myrank," opened file IELEP"
1990 call getenv_loc('SIDEPAR',sidename)
1991 open (isidep,file=sidename,status='old',action='read')
1992 c print *,"Processor",myrank," opened file ISIDEP"
1993 c print *,"Processor",myrank," opened parameter files"
1995 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1996 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1997 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1998 C Get parameter filenames and open the parameter files.
1999 call getenv_loc('BONDPAR',bondname)
2000 open (ibond,file=bondname,status='old')
2001 call getenv_loc('THETPAR',thetname)
2002 open (ithep,file=thetname,status='old')
2003 call getenv_loc('ROTPAR',rotname)
2004 open (irotam,file=rotname,status='old')
2005 call getenv_loc('TORPAR',torname)
2006 open (itorp,file=torname,status='old')
2007 call getenv_loc('TORDPAR',tordname)
2008 open (itordp,file=tordname,status='old')
2009 call getenv_loc('SCCORPAR',sccorname)
2010 open (isccor,file=sccorname,status='old')
2011 call getenv_loc('FOURIER',fouriername)
2012 open (ifourier,file=fouriername,status='old')
2013 call getenv_loc('ELEPAR',elename)
2014 open (ielep,file=elename,status='old')
2015 call getenv_loc('SIDEPAR',sidename)
2016 open (isidep,file=sidename,status='old')
2018 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2020 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2021 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2022 C Get parameter filenames and open the parameter files.
2023 call getenv_loc('BONDPAR',bondname)
2024 open (ibond,file=bondname,status='old',readonly)
2025 call getenv_loc('THETPAR',thetname)
2026 open (ithep,file=thetname,status='old',readonly)
2027 call getenv_loc('ROTPAR',rotname)
2028 open (irotam,file=rotname,status='old',readonly)
2029 call getenv_loc('TORPAR',torname)
2030 open (itorp,file=torname,status='old',readonly)
2031 call getenv_loc('TORDPAR',tordname)
2032 open (itordp,file=tordname,status='old',readonly)
2033 call getenv_loc('SCCORPAR',sccorname)
2034 open (isccor,file=sccorname,status='old',readonly)
2036 call getenv_loc('THETPARPDB',thetname_pdb)
2037 print *,"thetname_pdb ",thetname_pdb
2038 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2039 print *,ithep_pdb," opened"
2041 call getenv_loc('FOURIER',fouriername)
2042 open (ifourier,file=fouriername,status='old',readonly)
2043 call getenv_loc('ELEPAR',elename)
2044 open (ielep,file=elename,status='old',readonly)
2045 call getenv_loc('SIDEPAR',sidename)
2046 open (isidep,file=sidename,status='old',readonly)
2048 call getenv_loc('ROTPARPDB',rotname_pdb)
2049 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2054 C 8/9/01 In the newest version SCp interaction constants are read from a file
2055 C Use -DOLDSCP to use hard-coded constants instead.
2057 call getenv_loc('SCPPAR',scpname)
2058 #if defined(WINIFL) || defined(WINPGI)
2059 open (iscpp,file=scpname,status='old',readonly,shared)
2060 #elif (defined CRAY) || (defined AIX)
2061 open (iscpp,file=scpname,status='old',action='read')
2063 open (iscpp,file=scpname,status='old')
2065 open (iscpp,file=scpname,status='old',readonly)
2068 call getenv_loc('PATTERN',patname)
2069 #if defined(WINIFL) || defined(WINPGI)
2070 open (icbase,file=patname,status='old',readonly,shared)
2071 #elif (defined CRAY) || (defined AIX)
2072 open (icbase,file=patname,status='old',action='read')
2074 open (icbase,file=patname,status='old')
2076 open (icbase,file=patname,status='old',readonly)
2079 C Open output file only for CG processes
2080 c print *,"Processor",myrank," fg_rank",fg_rank
2081 if (fg_rank.eq.0) then
2083 if (nodes.eq.1) then
2086 npos = dlog10(dfloat(nodes-1))+1
2088 if (npos.lt.3) npos=3
2089 write (liczba,'(i1)') npos
2090 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2092 write (liczba,form) me
2093 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2094 & liczba(:ilen(liczba))
2095 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2097 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2099 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2100 & liczba(:ilen(liczba))//'.mol2'
2101 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2102 & liczba(:ilen(liczba))//'.stat'
2104 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2105 & //liczba(:ilen(liczba))//'.stat')
2106 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2109 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2110 & liczba(:ilen(liczba))//'.const'
2115 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2116 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2117 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2118 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2119 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2121 & call copy_to_tmp(pref_orig(:ile(pref_orig))//'_'//pot(:lenpot)//
2123 rest2name=prefix(:ilen(prefix))//'.rst'
2125 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2128 #if defined(AIX) || defined(PGI)
2129 if (me.eq.king .or. .not. out1file)
2130 & open(iout,file=outname,status='unknown')
2132 if (fg_rank.gt.0) then
2133 write (liczba,'(i3.3)') myrank/nfgtasks
2134 write (ll,'(bz,i3.3)') fg_rank
2135 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2140 open(igeom,file=intname,status='unknown',position='append')
2141 open(ipdb,file=pdbname,status='unknown')
2142 open(imol2,file=mol2name,status='unknown')
2143 open(istat,file=statname,status='unknown',position='append')
2145 c1out open(iout,file=outname,status='unknown')
2148 if (me.eq.king .or. .not.out1file)
2149 & open(iout,file=outname,status='unknown')
2151 if (fg_rank.gt.0) then
2152 write (liczba,'(i3.3)') myrank/nfgtasks
2153 write (ll,'(bz,i3.3)') fg_rank
2154 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2159 open(igeom,file=intname,status='unknown',access='append')
2160 open(ipdb,file=pdbname,status='unknown')
2161 open(imol2,file=mol2name,status='unknown')
2162 open(istat,file=statname,status='unknown',access='append')
2164 c1out open(iout,file=outname,status='unknown')
2167 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2168 csa_seed=prefix(:lenpre)//'.CSA.seed'
2169 csa_history=prefix(:lenpre)//'.CSA.history'
2170 csa_bank=prefix(:lenpre)//'.CSA.bank'
2171 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2172 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2173 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2174 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2175 csa_int=prefix(:lenpre)//'.int'
2176 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2177 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2178 csa_in=prefix(:lenpre)//'.CSA.in'
2179 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2182 write (iout,'(80(1h-))')
2183 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2184 write (iout,'(80(1h-))')
2185 write (iout,*) "Input file : ",
2186 & pref_orig(:ilen(pref_orig))//'.inp'
2187 write (iout,*) "Output file : ",
2188 & outname(:ilen(outname))
2190 write (iout,*) "Sidechain potential file : ",
2191 & sidename(:ilen(sidename))
2193 write (iout,*) "SCp potential file : ",
2194 & scpname(:ilen(scpname))
2196 write (iout,*) "Electrostatic potential file : ",
2197 & elename(:ilen(elename))
2198 write (iout,*) "Cumulant coefficient file : ",
2199 & fouriername(:ilen(fouriername))
2200 write (iout,*) "Torsional parameter file : ",
2201 & torname(:ilen(torname))
2202 write (iout,*) "Double torsional parameter file : ",
2203 & tordname(:ilen(tordname))
2204 write (iout,*) "SCCOR parameter file : ",
2205 & sccorname(:ilen(sccorname))
2206 write (iout,*) "Bond & inertia constant file : ",
2207 & bondname(:ilen(bondname))
2208 write (iout,*) "Bending parameter file : ",
2209 & thetname(:ilen(thetname))
2210 write (iout,*) "Rotamer parameter file : ",
2211 & rotname(:ilen(rotname))
2212 write (iout,*) "Thetpdb parameter file : ",
2213 & thetname_pdb(:ilen(thetname_pdb))
2214 write (iout,*) "Threading database : ",
2215 & patname(:ilen(patname))
2217 &write (iout,*)" DIRTMP : ",
2219 write (iout,'(80(1h-))')
2223 c----------------------------------------------------------------------------
2224 subroutine card_concat(card)
2225 implicit real*8 (a-h,o-z)
2226 include 'DIMENSIONS'
2227 include 'COMMON.IOUNITS'
2229 character*80 karta,ucase
2231 read (inp,'(a)') karta
2234 do while (karta(80:80).eq.'&')
2235 card=card(:ilen(card)+1)//karta(:79)
2236 read (inp,'(a)') karta
2239 card=card(:ilen(card)+1)//karta
2242 c----------------------------------------------------------------------------------
2244 implicit real*8 (a-h,o-z)
2245 include 'DIMENSIONS'
2246 include 'COMMON.CHAIN'
2247 include 'COMMON.IOUNITS'
2249 open(irest2,file=rest2name,status='unknown')
2250 read(irest2,*) totT,EK,potE,totE,t_bath
2252 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2255 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2258 read (irest2,*) iset
2263 c---------------------------------------------------------------------------------
2264 subroutine read_fragments
2265 implicit real*8 (a-h,o-z)
2266 include 'DIMENSIONS'
2270 include 'COMMON.SETUP'
2271 include 'COMMON.CHAIN'
2272 include 'COMMON.IOUNITS'
2274 include 'COMMON.CONTROL'
2275 read(inp,*) nset,nfrag,npair,nfrag_back
2276 if(me.eq.king.or..not.out1file)
2277 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2278 & " nfrag_back",nfrag_back
2280 read(inp,*) mset(iset)
2282 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2284 if(me.eq.king.or..not.out1file)
2285 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2286 & ifrag(2,i,iset), qinfrag(i,iset)
2289 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2291 if(me.eq.king.or..not.out1file)
2292 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2293 & ipair(2,i,iset), qinpair(i,iset)
2296 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2297 & wfrag_back(3,i,iset),
2298 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2299 if(me.eq.king.or..not.out1file)
2300 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2301 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2306 c-------------------------------------------------------------------------------
2307 subroutine read_dist_constr
2308 implicit real*8 (a-h,o-z)
2309 include 'DIMENSIONS'
2313 include 'COMMON.SETUP'
2314 include 'COMMON.CONTROL'
2315 include 'COMMON.CHAIN'
2316 include 'COMMON.IOUNITS'
2317 include 'COMMON.SBRIDGE'
2318 integer ifrag_(2,100),ipair_(2,100)
2319 double precision wfrag_(100),wpair_(100)
2320 character*500 controlcard
2322 write (iout,*) "Calling read_dist_constr"
2323 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2325 call card_concat(controlcard)
2326 call readi(controlcard,"NFRAG",nfrag_,0)
2327 call readi(controlcard,"NPAIR",npair_,0)
2328 call readi(controlcard,"NDIST",ndist_,0)
2329 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2330 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2331 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2332 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2333 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2334 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2335 c write (iout,*) "IFRAG"
2337 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2339 c write (iout,*) "IPAIR"
2341 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2345 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2346 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2347 & ifrag_(2,i)=nstart_sup+nsup-1
2348 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2350 if (wfrag_(i).gt.0.0d0) then
2351 do j=ifrag_(1,i),ifrag_(2,i)-1
2352 do k=j+1,ifrag_(2,i)
2353 c write (iout,*) "j",j," k",k
2355 if (constr_dist.eq.1) then
2360 forcon(nhpb)=wfrag_(i)
2361 else if (constr_dist.eq.2) then
2362 if (ddjk.le.dist_cut) then
2367 forcon(nhpb)=wfrag_(i)
2374 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2377 if (.not.out1file .or. me.eq.king)
2378 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2379 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2381 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2382 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2389 if (wpair_(i).gt.0.0d0) then
2397 do j=ifrag_(1,ii),ifrag_(2,ii)
2398 do k=ifrag_(1,jj),ifrag_(2,jj)
2402 forcon(nhpb)=wpair_(i)
2403 dhpb(nhpb)=dist(j,k)
2405 if (.not.out1file .or. me.eq.king)
2406 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2407 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2409 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2410 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2418 if (constr_dist.eq.11) then
2419 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2420 & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
2421 fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2424 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2425 & ibecarb(i),forcon(nhpb+1)
2427 if (forcon(nhpb+1).gt.0.0d0) then
2429 if (ibecarb(i).gt.0) then
2430 ihpb(i)=ihpb(i)+nres
2431 jhpb(i)=jhpb(i)+nres
2433 if (dhpb(nhpb).eq.0.0d0)
2434 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2436 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2437 C if (forcon(nhpb+1).gt.0.0d0) then
2439 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2441 if (.not.out1file .or. me.eq.king)
2442 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2443 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2445 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2446 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2453 c-------------------------------------------------------------------------------
2455 subroutine flush(iu)
2460 subroutine flush(iu)
2465 c------------------------------------------------------------------------------
2466 subroutine copy_to_tmp(source)
2467 include "DIMENSIONS"
2468 include "COMMON.IOUNITS"
2469 character*(*) source
2470 character* 256 tmpfile
2474 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2475 inquire(file=tmpfile,exist=ex)
2477 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2478 & " to temporary directory..."
2479 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2480 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2484 c------------------------------------------------------------------------------
2485 subroutine move_from_tmp(source)
2486 include "DIMENSIONS"
2487 include "COMMON.IOUNITS"
2488 character*(*) source
2491 write (*,*) "Moving ",source(:ilen(source)),
2492 & " from temporary directory to working directory"
2493 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2494 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2497 c------------------------------------------------------------------------------
2498 subroutine random_init(seed)
2500 C Initialize random number generator
2502 implicit real*8 (a-h,o-z)
2503 include 'DIMENSIONS'
2506 logical OKRandom, prng_restart
2508 integer iseed_array(4)
2510 include 'COMMON.IOUNITS'
2511 include 'COMMON.TIME1'
2512 include 'COMMON.THREAD'
2513 include 'COMMON.SBRIDGE'
2514 include 'COMMON.CONTROL'
2515 include 'COMMON.MCM'
2516 include 'COMMON.MAP'
2517 include 'COMMON.HEADER'
2518 include 'COMMON.CSA'
2519 include 'COMMON.CHAIN'
2520 include 'COMMON.MUCA'
2522 include 'COMMON.FFIELD'
2523 include 'COMMON.SETUP'
2524 iseed=-dint(dabs(seed))
2525 if (iseed.eq.0) then
2526 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2527 & 'Random seed undefined. The program will stop.'
2528 write (*,'(/80(1h*)/20x,a/80(1h*))')
2529 & 'Random seed undefined. The program will stop.'
2531 call mpi_finalize(mpi_comm_world,ierr)
2533 stop 'Bad random seed.'
2536 if (fg_rank.eq.0) then
2540 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2541 OKRandom = prng_restart(me,iseed)
2544 tmp=65536.0d0**(4-i)
2545 iseed_array(i) = dint(seed/tmp)
2546 seed=seed-iseed_array(i)*tmp
2549 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2550 & (iseed_array(i),i=1,4)
2551 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2552 & (iseed_array(i),i=1,4)
2553 OKRandom = prng_restart(me,iseed_array)
2556 c r1 = prng_next(me)
2557 r1=ran_number(0.0D0,1.0D0)
2559 & write (iout,*) 'ran_num',r1
2560 if (r1.lt.0.0d0) OKRandom=.false.
2562 if (.not.OKRandom) then
2563 write (iout,*) 'PRNG IS NOT WORKING!!!'
2564 print *,'PRNG IS NOT WORKING!!!'
2567 call mpi_abort(mpi_comm_world,error_msg,ierr)
2570 write (iout,*) 'too many processors for parallel prng'
2571 write (*,*) 'too many processors for parallel prng'
2579 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)