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 this can be applied
860 C for generate random conformation but is not implemented in this way
865 C begin reading theta constrains this is quartic constrains allowing to
866 C have smooth second derivative
867 if (with_theta_constr) then
868 C with_theta_constr is keyword allowing for occurance of theta constrains
869 read (inp,*) ntheta_constr
870 C ntheta_constr is the number of theta constrains
871 if (ntheta_constr.gt.0) then
873 read (inp,*) (itheta_constr(i),theta_constr0(i),
874 & theta_drange(i),for_thet_constr(i),
876 C the above code reads from 1 to ntheta_constr
877 C itheta_constr(i) residue i for which is theta_constr
878 C theta_constr0 the global minimum value
879 C theta_drange is range for which there is no energy penalty
880 C for_thet_constr is the force constant for quartic energy penalty
882 if(me.eq.king.or..not.out1file)then
884 & 'There are',ntheta_constr,' constraints on phi angles.'
886 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
892 theta_constr0(i)=deg2rad*theta_constr0(i)
893 theta_drange(i)=deg2rad*theta_drange(i)
895 C if(me.eq.king.or..not.out1file)
896 C & write (iout,*) 'FTORS',ftors
897 C do i=1,ntheta_constr
898 C ii = itheta_constr(i)
899 C thetabound(1,ii) = phi0(i)-drange(i)
900 C thetabound(2,ii) = phi0(i)+drange(i)
902 endif ! ntheta_constr.gt.0
903 endif! with_theta_constr
905 C with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
906 C write (iout,*) "with_dihed_constr ",with_dihed_constr
911 write (iout,'(a)') 'Boundaries in phi angle sampling:'
913 write (iout,'(a3,i5,2f10.1)')
914 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
920 cd print *,'NNT=',NNT,' NCT=',NCT
921 if (itype(1).eq.ntyp1) nnt=2
922 if (itype(nres).eq.ntyp1) nct=nct-1
924 if(me.eq.king.or..not.out1file)
925 & write (iout,'(a,i3)') 'nsup=',nsup
927 if (nsup.le.(nct-nnt+1)) then
928 do i=0,nct-nnt+1-nsup
929 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
935 & 'Error - sequences to be superposed do not match.'
938 do i=0,nsup-(nct-nnt+1)
939 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
941 nstart_sup=nstart_sup+i
947 & 'Error - sequences to be superposed do not match.'
950 if (nsup.eq.0) nsup=nct-nnt
951 if (nstart_sup.eq.0) nstart_sup=nnt
952 if (nstart_seq.eq.0) nstart_seq=nnt
953 if(me.eq.king.or..not.out1file)
954 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
955 & ' nstart_seq=',nstart_seq
957 c--- Zscore rms -------
958 if (nz_start.eq.0) nz_start=nnt
959 if (nz_end.eq.0 .and. nsup.gt.0) then
961 else if (nz_end.eq.0) then
964 if(me.eq.king.or..not.out1file)then
965 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
966 write (iout,*) 'IZ_SC=',iz_sc
968 c----------------------
971 if (.not.pdbref) then
972 call read_angles(inp,*38)
974 38 write (iout,'(a)') 'Error reading reference structure.'
976 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
977 stop 'Error reading reference structure'
981 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
991 call contact(.true.,ncont_ref,icont_ref,co)
995 write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
997 if (constr_dist.gt.0) call read_dist_constr
998 write (iout,*) "After read_dist_constr nhpb",nhpb
1000 if(me.eq.king.or..not.out1file)
1001 & write (iout,*) 'Contact order:',co
1003 if(me.eq.king.or..not.out1file)
1004 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1007 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1009 if(me.eq.king.or..not.out1file)
1010 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1011 & icont_ref(1,i),' ',
1012 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1016 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1017 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1018 & modecalc.ne.10) then
1019 C If input structure hasn't been supplied from the PDB file read or generate
1021 if (iranconf.eq.0 .and. .not. extconf) then
1022 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1023 & write (iout,'(a)') 'Initial geometry will be read in.'
1025 read(inp,'(8f10.5)',end=36,err=36)
1026 & ((c(l,k),l=1,3),k=1,nres),
1027 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1028 write (iout,*) "Exit READ_CART"
1029 write (iout,'(8f10.5)')
1030 & ((c(l,k),l=1,3),k=1,nres),
1031 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1032 call int_from_cart1(.true.)
1033 write (iout,*) "Finish INT_TO_CART"
1036 dc(j,i)=c(j,i+1)-c(j,i)
1037 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1041 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1043 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1044 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1050 call read_angles(inp,*36)
1053 36 write (iout,'(a)') 'Error reading angle file.'
1055 call mpi_finalize( MPI_COMM_WORLD,IERR )
1057 stop 'Error reading angle file.'
1059 else if (extconf) then
1060 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1061 & write (iout,'(a)') 'Extended chain initial geometry.'
1063 theta(i)=90d0*deg2rad
1066 phi(i)=180d0*deg2rad
1069 alph(i)=110d0*deg2rad
1072 omeg(i)=-120d0*deg2rad
1073 if (itype(i).le.0) omeg(i)=-omeg(i)
1076 if(me.eq.king.or..not.out1file)
1077 & write (iout,'(a)') 'Random-generated initial geometry.'
1081 if (me.eq.king .or. fg_rank.eq.0 .and. (
1082 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1086 call gen_rand_conf(itmp,*30)
1088 30 write (iout,*) 'Failed to generate random conformation',
1089 & ', itrial=',itrial
1090 write (*,*) 'Processor:',me,
1091 & ' Failed to generate random conformation',
1101 write (iout,'(a,i3,a)') 'Processor:',me,
1102 & ' error in generating random conformation.'
1103 write (*,'(a,i3,a)') 'Processor:',me,
1104 & ' error in generating random conformation.'
1107 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1113 call gen_rand_conf(itmp,*30)
1115 30 write (iout,*) 'Failed to generate random conformation',
1116 & ', itrial=',itrial
1117 write (*,*) 'Failed to generate random conformation',
1118 & ', itrial=',itrial
1120 write (iout,'(a,i3,a)') 'Processor:',me,
1121 & ' error in generating random conformation.'
1122 write (*,'(a,i3,a)') 'Processor:',me,
1123 & ' error in generating random conformation.'
1128 elseif (modecalc.eq.4) then
1129 read (inp,'(a)') intinname
1130 open (intin,file=intinname,status='old',err=333)
1131 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1132 & write (iout,'(a)') 'intinname',intinname
1133 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1135 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1137 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1139 stop 'Error opening angle file.'
1143 C Generate distance constraints, if the PDB structure is to be regularized.
1144 if (nthread.gt.0) then
1145 call read_threadbase
1148 if (me.eq.king .or. .not. out1file)
1150 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1151 write (iout,'(/a,i3,a)')
1152 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1153 write (iout,'(20i4)') (iss(i),i=1,ns)
1155 write(iout,*)"Running with dynamic disulfide-bond formation"
1157 write (iout,'(/a/)') 'Pre-formed links are:'
1163 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1164 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1170 if (ns.gt.0.and.dyn_ss) then
1174 forcon(i-nss)=forcon(i)
1181 dyn_ss_mask(iss(i))=.true.
1184 if (i2ndstr.gt.0) call secstrp2dihc
1185 c call geom_to_var(nvar,x)
1186 c call etotal(energia(0))
1187 c call enerprint(energia(0))
1188 c call briefout(0,etot)
1190 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1191 cd write (iout,'(a)') 'Variable list:'
1192 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1194 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1195 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1196 & 'Processor',myrank,': end reading molecular data.'
1201 c--------------------------------------------------------------------------
1202 logical function seq_comp(itypea,itypeb,length)
1204 integer length,itypea(length),itypeb(length)
1207 if (itypea(i).ne.itypeb(i)) then
1215 c-----------------------------------------------------------------------------
1216 subroutine read_bridge
1217 C Read information about disulfide bridges.
1218 implicit real*8 (a-h,o-z)
1219 include 'DIMENSIONS'
1223 include 'COMMON.IOUNITS'
1224 include 'COMMON.GEO'
1225 include 'COMMON.VAR'
1226 include 'COMMON.INTERACT'
1227 include 'COMMON.LOCAL'
1228 include 'COMMON.NAMES'
1229 include 'COMMON.CHAIN'
1230 include 'COMMON.FFIELD'
1231 include 'COMMON.SBRIDGE'
1232 include 'COMMON.HEADER'
1233 include 'COMMON.CONTROL'
1234 include 'COMMON.DBASE'
1235 include 'COMMON.THREAD'
1236 include 'COMMON.TIME1'
1237 include 'COMMON.SETUP'
1238 C Read bridging residues.
1239 read (inp,*) ns,(iss(i),i=1,ns)
1241 if(me.eq.king.or..not.out1file)
1242 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1243 C Check whether the specified bridging residues are cystines.
1245 if (itype(iss(i)).ne.1) then
1246 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1247 & 'Do you REALLY think that the residue ',
1248 & restyp(itype(iss(i))),i,
1249 & ' can form a disulfide bridge?!!!'
1250 write (*,'(2a,i3,a)')
1251 & 'Do you REALLY think that the residue ',
1252 & restyp(itype(iss(i))),i,
1253 & ' can form a disulfide bridge?!!!'
1255 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1260 C Read preformed bridges.
1262 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1264 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1267 C Check if the residues involved in bridges are in the specified list of
1268 C bridging residues.
1271 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1272 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1273 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1274 & ' contains residues present in other pairs.'
1275 write (*,'(a,i3,a)') 'Disulfide pair',i,
1276 & ' contains residues present in other pairs.'
1278 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1284 if (ihpb(i).eq.iss(j)) goto 10
1286 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1289 if (jhpb(i).eq.iss(j)) goto 20
1291 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1297 ihpb(i)=ihpb(i)+nres
1298 jhpb(i)=jhpb(i)+nres
1304 c----------------------------------------------------------------------------
1305 subroutine read_x(kanal,*)
1306 implicit real*8 (a-h,o-z)
1307 include 'DIMENSIONS'
1308 include 'COMMON.GEO'
1309 include 'COMMON.VAR'
1310 include 'COMMON.CHAIN'
1311 include 'COMMON.IOUNITS'
1312 include 'COMMON.CONTROL'
1313 include 'COMMON.LOCAL'
1314 include 'COMMON.INTERACT'
1315 c Read coordinates from input
1317 read(kanal,'(8f10.5)',end=10,err=10)
1318 & ((c(l,k),l=1,3),k=1,nres),
1319 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1322 c(j,2*nres)=c(j,nres)
1324 call int_from_cart1(.false.)
1327 dc(j,i)=c(j,i+1)-c(j,i)
1328 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1332 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1334 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1335 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1343 c----------------------------------------------------------------------------
1344 subroutine read_threadbase
1345 implicit real*8 (a-h,o-z)
1346 include 'DIMENSIONS'
1347 include 'COMMON.IOUNITS'
1348 include 'COMMON.GEO'
1349 include 'COMMON.VAR'
1350 include 'COMMON.INTERACT'
1351 include 'COMMON.LOCAL'
1352 include 'COMMON.NAMES'
1353 include 'COMMON.CHAIN'
1354 include 'COMMON.FFIELD'
1355 include 'COMMON.SBRIDGE'
1356 include 'COMMON.HEADER'
1357 include 'COMMON.CONTROL'
1358 include 'COMMON.DBASE'
1359 include 'COMMON.THREAD'
1360 include 'COMMON.TIME1'
1361 C Read pattern database for threading.
1362 read (icbase,*) nseq
1364 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1365 & nres_base(2,i),nres_base(3,i)
1366 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1368 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1369 c & nres_base(2,i),nres_base(3,i)
1370 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1374 if (weidis.eq.0.0D0) weidis=0.1D0
1383 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1384 write (iout,'(a,i5)') 'nexcl: ',nexcl
1385 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1388 c------------------------------------------------------------------------------
1389 subroutine setup_var
1390 implicit real*8 (a-h,o-z)
1391 include 'DIMENSIONS'
1392 include 'COMMON.IOUNITS'
1393 include 'COMMON.GEO'
1394 include 'COMMON.VAR'
1395 include 'COMMON.INTERACT'
1396 include 'COMMON.LOCAL'
1397 include 'COMMON.NAMES'
1398 include 'COMMON.CHAIN'
1399 include 'COMMON.FFIELD'
1400 include 'COMMON.SBRIDGE'
1401 include 'COMMON.HEADER'
1402 include 'COMMON.CONTROL'
1403 include 'COMMON.DBASE'
1404 include 'COMMON.THREAD'
1405 include 'COMMON.TIME1'
1406 C Set up variable list.
1412 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1414 ialph(i,1)=nvar+nside
1418 if (indphi.gt.0) then
1420 else if (indback.gt.0) then
1425 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1428 c----------------------------------------------------------------------------
1429 subroutine gen_dist_constr
1430 C Generate CA distance constraints.
1431 implicit real*8 (a-h,o-z)
1432 include 'DIMENSIONS'
1433 include 'COMMON.IOUNITS'
1434 include 'COMMON.GEO'
1435 include 'COMMON.VAR'
1436 include 'COMMON.INTERACT'
1437 include 'COMMON.LOCAL'
1438 include 'COMMON.NAMES'
1439 include 'COMMON.CHAIN'
1440 include 'COMMON.FFIELD'
1441 include 'COMMON.SBRIDGE'
1442 include 'COMMON.HEADER'
1443 include 'COMMON.CONTROL'
1444 include 'COMMON.DBASE'
1445 include 'COMMON.THREAD'
1446 include 'COMMON.TIME1'
1447 dimension itype_pdb(maxres)
1448 common /pizda/ itype_pdb
1450 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1451 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1452 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1454 do i=nstart_sup,nstart_sup+nsup-1
1455 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1456 cd & ' seq_pdb', restyp(itype_pdb(i))
1457 do j=i+2,nstart_sup+nsup-1
1459 ihpb(nhpb)=i+nstart_seq-nstart_sup
1460 jhpb(nhpb)=j+nstart_seq-nstart_sup
1462 dhpb(nhpb)=dist(i,j)
1465 cd write (iout,'(a)') 'Distance constraints:'
1470 cd if (ii.gt.nres) then
1475 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1476 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1477 cd & dhpb(i),forcon(i)
1481 c----------------------------------------------------------------------------
1483 implicit real*8 (a-h,o-z)
1484 include 'DIMENSIONS'
1485 include 'COMMON.MAP'
1486 include 'COMMON.IOUNITS'
1487 character*3 angid(4) /'THE','PHI','ALP','OME'/
1488 character*80 mapcard,ucase
1490 read (inp,'(a)') mapcard
1491 mapcard=ucase(mapcard)
1492 if (index(mapcard,'PHI').gt.0) then
1494 else if (index(mapcard,'THE').gt.0) then
1496 else if (index(mapcard,'ALP').gt.0) then
1498 else if (index(mapcard,'OME').gt.0) then
1501 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1502 stop 'Error - illegal variable spec in MAP card.'
1504 call readi (mapcard,'RES1',res1(imap),0)
1505 call readi (mapcard,'RES2',res2(imap),0)
1506 if (res1(imap).eq.0) then
1507 res1(imap)=res2(imap)
1508 else if (res2(imap).eq.0) then
1509 res2(imap)=res1(imap)
1511 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1513 & 'Error - illegal definition of variable group in MAP.'
1514 stop 'Error - illegal definition of variable group in MAP.'
1516 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1517 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1518 call readi(mapcard,'NSTEP',nstep(imap),0)
1519 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1521 & 'Illegal boundary and/or step size specification in MAP.'
1522 stop 'Illegal boundary and/or step size specification in MAP.'
1527 c----------------------------------------------------------------------------
1529 implicit real*8 (a-h,o-z)
1530 include 'DIMENSIONS'
1531 include 'COMMON.IOUNITS'
1532 include 'COMMON.GEO'
1533 include 'COMMON.CSA'
1534 include 'COMMON.BANK'
1535 include 'COMMON.CONTROL'
1537 character*620 mcmcard
1538 call card_concat(mcmcard)
1540 call readi(mcmcard,'NCONF',nconf,50)
1541 call readi(mcmcard,'NADD',nadd,0)
1542 call readi(mcmcard,'JSTART',jstart,1)
1543 call readi(mcmcard,'JEND',jend,1)
1544 call readi(mcmcard,'NSTMAX',nstmax,500000)
1545 call readi(mcmcard,'N0',n0,1)
1546 call readi(mcmcard,'N1',n1,6)
1547 call readi(mcmcard,'N2',n2,4)
1548 call readi(mcmcard,'N3',n3,0)
1549 call readi(mcmcard,'N4',n4,0)
1550 call readi(mcmcard,'N5',n5,0)
1551 call readi(mcmcard,'N6',n6,10)
1552 call readi(mcmcard,'N7',n7,0)
1553 call readi(mcmcard,'N8',n8,0)
1554 call readi(mcmcard,'N9',n9,0)
1555 call readi(mcmcard,'N14',n14,0)
1556 call readi(mcmcard,'N15',n15,0)
1557 call readi(mcmcard,'N16',n16,0)
1558 call readi(mcmcard,'N17',n17,0)
1559 call readi(mcmcard,'N18',n18,0)
1561 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1563 call readi(mcmcard,'NDIFF',ndiff,2)
1564 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1565 call readi(mcmcard,'IS1',is1,1)
1566 call readi(mcmcard,'IS2',is2,8)
1567 call readi(mcmcard,'NRAN0',nran0,4)
1568 call readi(mcmcard,'NRAN1',nran1,2)
1569 call readi(mcmcard,'IRR',irr,1)
1570 call readi(mcmcard,'NSEED',nseed,20)
1571 call readi(mcmcard,'NTOTAL',ntotal,10000)
1572 call reada(mcmcard,'CUT1',cut1,2.0d0)
1573 call reada(mcmcard,'CUT2',cut2,5.0d0)
1574 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1575 call readi(mcmcard,'ICMAX',icmax,3)
1576 call readi(mcmcard,'IRESTART',irestart,0)
1577 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1580 call reada(mcmcard,'DELE',dele,20.0d0)
1581 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1582 call readi(mcmcard,'IREF',iref,0)
1583 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1584 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1585 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1586 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1587 write (iout,*) "NCONF_IN",nconf_in
1590 c----------------------------------------------------------------------------
1591 cfmc subroutine mcmfread
1592 cfmc implicit real*8 (a-h,o-z)
1593 cfmc include 'DIMENSIONS'
1594 cfmc include 'COMMON.MCMF'
1595 cfmc include 'COMMON.IOUNITS'
1596 cfmc include 'COMMON.GEO'
1597 cfmc character*80 ucase
1598 cfmc character*620 mcmcard
1599 cfmc call card_concat(mcmcard)
1601 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1602 cfmc write(iout,*)'MAXRANT=',maxrant
1603 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1604 cfmc write(iout,*)'MAXFAM=',maxfam
1605 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1606 cfmc write(iout,*)'NNET1=',nnet1
1607 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1608 cfmc write(iout,*)'NNET2=',nnet2
1609 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1610 cfmc write(iout,*)'NNET3=',nnet3
1611 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1612 cfmc write(iout,*)'ILASTT=',ilastt
1613 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1614 cfmc write(iout,*)'MAXSTR=',maxstr
1615 cfmc maxstr_f=maxstr/maxfam
1616 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1617 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1618 cfmc write(iout,*)'NMCMF=',nmcmf
1619 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1620 cfmc write(iout,*)'IFOCUS=',ifocus
1621 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1622 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1623 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1624 cfmc write(iout,*)'INTPRT=',intprt
1625 cfmc call readi(mcmcard,'IPRT',iprt,100)
1626 cfmc write(iout,*)'IPRT=',iprt
1627 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1628 cfmc write(iout,*)'IMAXTR=',imaxtr
1629 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1630 cfmc write(iout,*)'MAXEVEN=',maxeven
1631 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1632 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1633 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1634 cfmc write(iout,*)'INIMIN=',inimin
1635 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1636 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1637 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1638 cfmc write(iout,*)'NTHREAD=',nthread
1639 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1640 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1641 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1642 cfmc write(iout,*)'MAXPERT=',maxpert
1643 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1644 cfmc write(iout,*)'IRMSD=',irmsd
1645 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1646 cfmc write(iout,*)'DENEMIN=',denemin
1647 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1648 cfmc write(iout,*)'RCUT1S=',rcut1s
1649 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1650 cfmc write(iout,*)'RCUT1E=',rcut1e
1651 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1652 cfmc write(iout,*)'RCUT2S=',rcut2s
1653 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1654 cfmc write(iout,*)'RCUT2E=',rcut2e
1655 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1656 cfmc write(iout,*)'DPERT1=',d_pert1
1657 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1658 cfmc write(iout,*)'DPERT1A=',d_pert1a
1659 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1660 cfmc write(iout,*)'DPERT2=',d_pert2
1661 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1662 cfmc write(iout,*)'DPERT2A=',d_pert2a
1663 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1664 cfmc write(iout,*)'DPERT2B=',d_pert2b
1665 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1666 cfmc write(iout,*)'DPERT2C=',d_pert2c
1667 cfmc d_pert1=deg2rad*d_pert1
1668 cfmc d_pert1a=deg2rad*d_pert1a
1669 cfmc d_pert2=deg2rad*d_pert2
1670 cfmc d_pert2a=deg2rad*d_pert2a
1671 cfmc d_pert2b=deg2rad*d_pert2b
1672 cfmc d_pert2c=deg2rad*d_pert2c
1673 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1674 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1675 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1676 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1677 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1678 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1679 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1680 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1681 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1682 cfmc write(iout,*)'RCUTINI=',rcutini
1683 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1684 cfmc write(iout,*)'GRAT=',grat
1685 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1686 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1690 c----------------------------------------------------------------------------
1692 implicit real*8 (a-h,o-z)
1693 include 'DIMENSIONS'
1694 include 'COMMON.MCM'
1695 include 'COMMON.MCE'
1696 include 'COMMON.IOUNITS'
1698 character*320 mcmcard
1699 call card_concat(mcmcard)
1700 call readi(mcmcard,'MAXACC',maxacc,100)
1701 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1702 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1703 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1704 call readi(mcmcard,'MAXREPM',maxrepm,200)
1705 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1706 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1707 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1708 call reada(mcmcard,'E_UP',e_up,5.0D0)
1709 call reada(mcmcard,'DELTE',delte,0.1D0)
1710 call readi(mcmcard,'NSWEEP',nsweep,5)
1711 call readi(mcmcard,'NSTEPH',nsteph,0)
1712 call readi(mcmcard,'NSTEPC',nstepc,0)
1713 call reada(mcmcard,'TMIN',tmin,298.0D0)
1714 call reada(mcmcard,'TMAX',tmax,298.0D0)
1715 call readi(mcmcard,'NWINDOW',nwindow,0)
1716 call readi(mcmcard,'PRINT_MC',print_mc,0)
1717 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1718 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1719 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1720 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1721 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1722 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1723 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1724 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1725 if (nwindow.gt.0) then
1726 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1728 winlen(i)=winend(i)-winstart(i)+1
1731 if (tmax.lt.tmin) tmax=tmin
1732 if (tmax.eq.tmin) then
1736 if (nstepc.gt.0 .and. nsteph.gt.0) then
1737 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1738 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1740 C Probabilities of different move types
1741 sumpro_type(0)=0.0D0
1742 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1743 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1744 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1745 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1746 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1747 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1748 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1750 print *,'i',i,' sumprotype',sumpro_type(i)
1751 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1752 print *,'i',i,' sumprotype',sumpro_type(i)
1756 c----------------------------------------------------------------------------
1757 subroutine read_minim
1758 implicit real*8 (a-h,o-z)
1759 include 'DIMENSIONS'
1760 include 'COMMON.MINIM'
1761 include 'COMMON.IOUNITS'
1763 character*320 minimcard
1764 call card_concat(minimcard)
1765 call readi(minimcard,'MAXMIN',maxmin,2000)
1766 call readi(minimcard,'MAXFUN',maxfun,5000)
1767 call readi(minimcard,'MINMIN',minmin,maxmin)
1768 call readi(minimcard,'MINFUN',minfun,maxmin)
1769 call reada(minimcard,'TOLF',tolf,1.0D-2)
1770 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1771 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1772 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1773 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1774 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1775 & 'Options in energy minimization:'
1776 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1777 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1778 & 'MinMin:',MinMin,' MinFun:',MinFun,
1779 & ' TolF:',TolF,' RTolF:',RTolF
1782 c----------------------------------------------------------------------------
1783 subroutine read_angles(kanal,*)
1784 implicit real*8 (a-h,o-z)
1785 include 'DIMENSIONS'
1786 include 'COMMON.GEO'
1787 include 'COMMON.VAR'
1788 include 'COMMON.CHAIN'
1789 include 'COMMON.IOUNITS'
1790 include 'COMMON.CONTROL'
1791 c Read angles from input
1793 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1794 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1795 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1796 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1799 c 9/7/01 avoid 180 deg valence angle
1800 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1802 theta(i)=deg2rad*theta(i)
1803 phi(i)=deg2rad*phi(i)
1804 alph(i)=deg2rad*alph(i)
1805 omeg(i)=deg2rad*omeg(i)
1810 c----------------------------------------------------------------------------
1811 subroutine reada(rekord,lancuch,wartosc,default)
1813 character*(*) rekord,lancuch
1814 double precision wartosc,default
1817 iread=index(rekord,lancuch)
1818 if (iread.eq.0) then
1822 iread=iread+ilen(lancuch)+1
1823 read (rekord(iread:),*,err=10,end=10) wartosc
1828 c----------------------------------------------------------------------------
1829 subroutine readi(rekord,lancuch,wartosc,default)
1831 character*(*) rekord,lancuch
1832 integer wartosc,default
1835 iread=index(rekord,lancuch)
1836 if (iread.eq.0) then
1840 iread=iread+ilen(lancuch)+1
1841 read (rekord(iread:),*,err=10,end=10) wartosc
1846 c----------------------------------------------------------------------------
1847 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1850 integer tablica(dim),default
1851 character*(*) rekord,lancuch
1858 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1859 if (iread.eq.0) return
1860 iread=iread+ilen(lancuch)+1
1861 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1864 c----------------------------------------------------------------------------
1865 subroutine multreada(rekord,lancuch,tablica,dim,default)
1868 double precision tablica(dim),default
1869 character*(*) rekord,lancuch
1876 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1877 if (iread.eq.0) return
1878 iread=iread+ilen(lancuch)+1
1879 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1882 c----------------------------------------------------------------------------
1883 subroutine openunits
1884 implicit real*8 (a-h,o-z)
1885 include 'DIMENSIONS'
1888 character*16 form,nodename
1891 include 'COMMON.SETUP'
1892 include 'COMMON.IOUNITS'
1894 include 'COMMON.CONTROL'
1895 integer lenpre,lenpot,ilen,lentmp
1897 character*3 out1file_text,ucase
1900 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1901 call getenv_loc("PREFIX",prefix)
1903 call getenv_loc("POT",pot)
1904 call getenv_loc("DIRTMP",tmpdir)
1905 call getenv_loc("CURDIR",curdir)
1906 call getenv_loc("OUT1FILE",out1file_text)
1907 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1908 out1file_text=ucase(out1file_text)
1909 if (out1file_text(1:1).eq."Y") then
1912 out1file=fg_rank.gt.0
1917 if (lentmp.gt.0) then
1918 write (*,'(80(1h!))')
1919 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1920 write (*,'(80(1h!))')
1921 write (*,*)"All output files will be on node /tmp directory."
1923 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1924 if (me.eq.king) then
1925 write (*,*) "The master node is ",nodename
1926 else if (fg_rank.eq.0) then
1927 write (*,*) "I am the CG slave node ",nodename
1929 write (*,*) "I am the FG slave node ",nodename
1932 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1933 lenpre = lentmp+lenpre+1
1935 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1936 C Get the names and open the input files
1937 #if defined(WINIFL) || defined(WINPGI)
1938 open(1,file=pref_orig(:ilen(pref_orig))//
1939 & '.inp',status='old',readonly,shared)
1940 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1941 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1942 C Get parameter filenames and open the parameter files.
1943 call getenv_loc('BONDPAR',bondname)
1944 open (ibond,file=bondname,status='old',readonly,shared)
1945 call getenv_loc('THETPAR',thetname)
1946 open (ithep,file=thetname,status='old',readonly,shared)
1947 call getenv_loc('ROTPAR',rotname)
1948 open (irotam,file=rotname,status='old',readonly,shared)
1949 call getenv_loc('TORPAR',torname)
1950 open (itorp,file=torname,status='old',readonly,shared)
1951 call getenv_loc('TORDPAR',tordname)
1952 open (itordp,file=tordname,status='old',readonly,shared)
1953 call getenv_loc('FOURIER',fouriername)
1954 open (ifourier,file=fouriername,status='old',readonly,shared)
1955 call getenv_loc('ELEPAR',elename)
1956 open (ielep,file=elename,status='old',readonly,shared)
1957 call getenv_loc('SIDEPAR',sidename)
1958 open (isidep,file=sidename,status='old',readonly,shared)
1959 #elif (defined CRAY) || (defined AIX)
1960 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1962 c print *,"Processor",myrank," opened file 1"
1963 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1964 c print *,"Processor",myrank," opened file 9"
1965 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1966 C Get parameter filenames and open the parameter files.
1967 call getenv_loc('BONDPAR',bondname)
1968 open (ibond,file=bondname,status='old',action='read')
1969 c print *,"Processor",myrank," opened file IBOND"
1970 call getenv_loc('THETPAR',thetname)
1971 open (ithep,file=thetname,status='old',action='read')
1972 c print *,"Processor",myrank," opened file ITHEP"
1973 call getenv_loc('ROTPAR',rotname)
1974 open (irotam,file=rotname,status='old',action='read')
1975 c print *,"Processor",myrank," opened file IROTAM"
1976 call getenv_loc('TORPAR',torname)
1977 open (itorp,file=torname,status='old',action='read')
1978 c print *,"Processor",myrank," opened file ITORP"
1979 call getenv_loc('TORDPAR',tordname)
1980 open (itordp,file=tordname,status='old',action='read')
1981 c print *,"Processor",myrank," opened file ITORDP"
1982 call getenv_loc('SCCORPAR',sccorname)
1983 open (isccor,file=sccorname,status='old',action='read')
1984 c print *,"Processor",myrank," opened file ISCCOR"
1985 call getenv_loc('FOURIER',fouriername)
1986 open (ifourier,file=fouriername,status='old',action='read')
1987 c print *,"Processor",myrank," opened file IFOURIER"
1988 call getenv_loc('ELEPAR',elename)
1989 open (ielep,file=elename,status='old',action='read')
1990 c print *,"Processor",myrank," opened file IELEP"
1991 call getenv_loc('SIDEPAR',sidename)
1992 open (isidep,file=sidename,status='old',action='read')
1993 c print *,"Processor",myrank," opened file ISIDEP"
1994 c print *,"Processor",myrank," opened parameter files"
1996 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1997 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1998 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1999 C Get parameter filenames and open the parameter files.
2000 call getenv_loc('BONDPAR',bondname)
2001 open (ibond,file=bondname,status='old')
2002 call getenv_loc('THETPAR',thetname)
2003 open (ithep,file=thetname,status='old')
2004 call getenv_loc('ROTPAR',rotname)
2005 open (irotam,file=rotname,status='old')
2006 call getenv_loc('TORPAR',torname)
2007 open (itorp,file=torname,status='old')
2008 call getenv_loc('TORDPAR',tordname)
2009 open (itordp,file=tordname,status='old')
2010 call getenv_loc('SCCORPAR',sccorname)
2011 open (isccor,file=sccorname,status='old')
2012 call getenv_loc('FOURIER',fouriername)
2013 open (ifourier,file=fouriername,status='old')
2014 call getenv_loc('ELEPAR',elename)
2015 open (ielep,file=elename,status='old')
2016 call getenv_loc('SIDEPAR',sidename)
2017 open (isidep,file=sidename,status='old')
2019 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2021 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2022 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2023 C Get parameter filenames and open the parameter files.
2024 call getenv_loc('BONDPAR',bondname)
2025 open (ibond,file=bondname,status='old',readonly)
2026 call getenv_loc('THETPAR',thetname)
2027 open (ithep,file=thetname,status='old',readonly)
2028 call getenv_loc('ROTPAR',rotname)
2029 open (irotam,file=rotname,status='old',readonly)
2030 call getenv_loc('TORPAR',torname)
2031 open (itorp,file=torname,status='old',readonly)
2032 call getenv_loc('TORDPAR',tordname)
2033 open (itordp,file=tordname,status='old',readonly)
2034 call getenv_loc('SCCORPAR',sccorname)
2035 open (isccor,file=sccorname,status='old',readonly)
2037 call getenv_loc('THETPARPDB',thetname_pdb)
2038 print *,"thetname_pdb ",thetname_pdb
2039 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2040 print *,ithep_pdb," opened"
2042 call getenv_loc('FOURIER',fouriername)
2043 open (ifourier,file=fouriername,status='old',readonly)
2044 call getenv_loc('ELEPAR',elename)
2045 open (ielep,file=elename,status='old',readonly)
2046 call getenv_loc('SIDEPAR',sidename)
2047 open (isidep,file=sidename,status='old',readonly)
2049 call getenv_loc('ROTPARPDB',rotname_pdb)
2050 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2055 C 8/9/01 In the newest version SCp interaction constants are read from a file
2056 C Use -DOLDSCP to use hard-coded constants instead.
2058 call getenv_loc('SCPPAR',scpname)
2059 #if defined(WINIFL) || defined(WINPGI)
2060 open (iscpp,file=scpname,status='old',readonly,shared)
2061 #elif (defined CRAY) || (defined AIX)
2062 open (iscpp,file=scpname,status='old',action='read')
2064 open (iscpp,file=scpname,status='old')
2066 open (iscpp,file=scpname,status='old',readonly)
2069 call getenv_loc('PATTERN',patname)
2070 #if defined(WINIFL) || defined(WINPGI)
2071 open (icbase,file=patname,status='old',readonly,shared)
2072 #elif (defined CRAY) || (defined AIX)
2073 open (icbase,file=patname,status='old',action='read')
2075 open (icbase,file=patname,status='old')
2077 open (icbase,file=patname,status='old',readonly)
2080 C Open output file only for CG processes
2081 c print *,"Processor",myrank," fg_rank",fg_rank
2082 if (fg_rank.eq.0) then
2084 if (nodes.eq.1) then
2087 npos = dlog10(dfloat(nodes-1))+1
2089 if (npos.lt.3) npos=3
2090 write (liczba,'(i1)') npos
2091 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2093 write (liczba,form) me
2094 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2095 & liczba(:ilen(liczba))
2096 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2098 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2100 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2101 & liczba(:ilen(liczba))//'.mol2'
2102 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2103 & liczba(:ilen(liczba))//'.stat'
2105 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2106 & //liczba(:ilen(liczba))//'.stat')
2107 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2110 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2111 & liczba(:ilen(liczba))//'.const'
2116 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2117 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2118 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2119 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2120 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2122 & call copy_to_tmp(pref_orig(:ile(pref_orig))//'_'//pot(:lenpot)//
2124 rest2name=prefix(:ilen(prefix))//'.rst'
2126 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2129 #if defined(AIX) || defined(PGI)
2130 if (me.eq.king .or. .not. out1file)
2131 & open(iout,file=outname,status='unknown')
2133 if (fg_rank.gt.0) then
2134 write (liczba,'(i3.3)') myrank/nfgtasks
2135 write (ll,'(bz,i3.3)') fg_rank
2136 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2141 open(igeom,file=intname,status='unknown',position='append')
2142 open(ipdb,file=pdbname,status='unknown')
2143 open(imol2,file=mol2name,status='unknown')
2144 open(istat,file=statname,status='unknown',position='append')
2146 c1out open(iout,file=outname,status='unknown')
2149 if (me.eq.king .or. .not.out1file)
2150 & open(iout,file=outname,status='unknown')
2152 if (fg_rank.gt.0) then
2153 write (liczba,'(i3.3)') myrank/nfgtasks
2154 write (ll,'(bz,i3.3)') fg_rank
2155 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2160 open(igeom,file=intname,status='unknown',access='append')
2161 open(ipdb,file=pdbname,status='unknown')
2162 open(imol2,file=mol2name,status='unknown')
2163 open(istat,file=statname,status='unknown',access='append')
2165 c1out open(iout,file=outname,status='unknown')
2168 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2169 csa_seed=prefix(:lenpre)//'.CSA.seed'
2170 csa_history=prefix(:lenpre)//'.CSA.history'
2171 csa_bank=prefix(:lenpre)//'.CSA.bank'
2172 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2173 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2174 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2175 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2176 csa_int=prefix(:lenpre)//'.int'
2177 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2178 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2179 csa_in=prefix(:lenpre)//'.CSA.in'
2180 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2183 write (iout,'(80(1h-))')
2184 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2185 write (iout,'(80(1h-))')
2186 write (iout,*) "Input file : ",
2187 & pref_orig(:ilen(pref_orig))//'.inp'
2188 write (iout,*) "Output file : ",
2189 & outname(:ilen(outname))
2191 write (iout,*) "Sidechain potential file : ",
2192 & sidename(:ilen(sidename))
2194 write (iout,*) "SCp potential file : ",
2195 & scpname(:ilen(scpname))
2197 write (iout,*) "Electrostatic potential file : ",
2198 & elename(:ilen(elename))
2199 write (iout,*) "Cumulant coefficient file : ",
2200 & fouriername(:ilen(fouriername))
2201 write (iout,*) "Torsional parameter file : ",
2202 & torname(:ilen(torname))
2203 write (iout,*) "Double torsional parameter file : ",
2204 & tordname(:ilen(tordname))
2205 write (iout,*) "SCCOR parameter file : ",
2206 & sccorname(:ilen(sccorname))
2207 write (iout,*) "Bond & inertia constant file : ",
2208 & bondname(:ilen(bondname))
2209 write (iout,*) "Bending parameter file : ",
2210 & thetname(:ilen(thetname))
2211 write (iout,*) "Rotamer parameter file : ",
2212 & rotname(:ilen(rotname))
2213 write (iout,*) "Thetpdb parameter file : ",
2214 & thetname_pdb(:ilen(thetname_pdb))
2215 write (iout,*) "Threading database : ",
2216 & patname(:ilen(patname))
2218 &write (iout,*)" DIRTMP : ",
2220 write (iout,'(80(1h-))')
2224 c----------------------------------------------------------------------------
2225 subroutine card_concat(card)
2226 implicit real*8 (a-h,o-z)
2227 include 'DIMENSIONS'
2228 include 'COMMON.IOUNITS'
2230 character*80 karta,ucase
2232 read (inp,'(a)') karta
2235 do while (karta(80:80).eq.'&')
2236 card=card(:ilen(card)+1)//karta(:79)
2237 read (inp,'(a)') karta
2240 card=card(:ilen(card)+1)//karta
2243 c----------------------------------------------------------------------------------
2245 implicit real*8 (a-h,o-z)
2246 include 'DIMENSIONS'
2247 include 'COMMON.CHAIN'
2248 include 'COMMON.IOUNITS'
2250 open(irest2,file=rest2name,status='unknown')
2251 read(irest2,*) totT,EK,potE,totE,t_bath
2253 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2256 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2259 read (irest2,*) iset
2264 c---------------------------------------------------------------------------------
2265 subroutine read_fragments
2266 implicit real*8 (a-h,o-z)
2267 include 'DIMENSIONS'
2271 include 'COMMON.SETUP'
2272 include 'COMMON.CHAIN'
2273 include 'COMMON.IOUNITS'
2275 include 'COMMON.CONTROL'
2276 read(inp,*) nset,nfrag,npair,nfrag_back
2277 if(me.eq.king.or..not.out1file)
2278 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2279 & " nfrag_back",nfrag_back
2281 read(inp,*) mset(iset)
2283 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2285 if(me.eq.king.or..not.out1file)
2286 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2287 & ifrag(2,i,iset), qinfrag(i,iset)
2290 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2292 if(me.eq.king.or..not.out1file)
2293 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2294 & ipair(2,i,iset), qinpair(i,iset)
2297 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2298 & wfrag_back(3,i,iset),
2299 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2300 if(me.eq.king.or..not.out1file)
2301 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2302 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2307 c-------------------------------------------------------------------------------
2308 subroutine read_dist_constr
2309 implicit real*8 (a-h,o-z)
2310 include 'DIMENSIONS'
2314 include 'COMMON.SETUP'
2315 include 'COMMON.CONTROL'
2316 include 'COMMON.CHAIN'
2317 include 'COMMON.IOUNITS'
2318 include 'COMMON.SBRIDGE'
2319 integer ifrag_(2,100),ipair_(2,100)
2320 double precision wfrag_(100),wpair_(100)
2321 character*500 controlcard
2323 write (iout,*) "Calling read_dist_constr"
2324 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2326 call card_concat(controlcard)
2327 call readi(controlcard,"NFRAG",nfrag_,0)
2328 call readi(controlcard,"NPAIR",npair_,0)
2329 call readi(controlcard,"NDIST",ndist_,0)
2330 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2331 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2332 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2333 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2334 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2335 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2336 c write (iout,*) "IFRAG"
2338 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2340 c write (iout,*) "IPAIR"
2342 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2346 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2347 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2348 & ifrag_(2,i)=nstart_sup+nsup-1
2349 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2351 if (wfrag_(i).gt.0.0d0) then
2352 do j=ifrag_(1,i),ifrag_(2,i)-1
2353 do k=j+1,ifrag_(2,i)
2354 c write (iout,*) "j",j," k",k
2356 if (constr_dist.eq.1) then
2361 forcon(nhpb)=wfrag_(i)
2362 else if (constr_dist.eq.2) then
2363 if (ddjk.le.dist_cut) then
2368 forcon(nhpb)=wfrag_(i)
2375 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2378 if (.not.out1file .or. me.eq.king)
2379 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2380 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2382 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2383 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2390 if (wpair_(i).gt.0.0d0) then
2398 do j=ifrag_(1,ii),ifrag_(2,ii)
2399 do k=ifrag_(1,jj),ifrag_(2,jj)
2403 forcon(nhpb)=wpair_(i)
2404 dhpb(nhpb)=dist(j,k)
2406 if (.not.out1file .or. me.eq.king)
2407 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2408 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2410 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2411 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2419 if (constr_dist.eq.11) then
2420 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2421 & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
2422 fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2425 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2426 & ibecarb(i),forcon(nhpb+1)
2428 if (forcon(nhpb+1).gt.0.0d0) then
2430 if (ibecarb(i).gt.0) then
2431 ihpb(i)=ihpb(i)+nres
2432 jhpb(i)=jhpb(i)+nres
2434 if (dhpb(nhpb).eq.0.0d0)
2435 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2437 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2438 C if (forcon(nhpb+1).gt.0.0d0) then
2440 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2442 if (.not.out1file .or. me.eq.king)
2443 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2444 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2446 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2447 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2454 c-------------------------------------------------------------------------------
2456 subroutine flush(iu)
2461 subroutine flush(iu)
2466 c------------------------------------------------------------------------------
2467 subroutine copy_to_tmp(source)
2468 include "DIMENSIONS"
2469 include "COMMON.IOUNITS"
2470 character*(*) source
2471 character* 256 tmpfile
2475 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2476 inquire(file=tmpfile,exist=ex)
2478 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2479 & " to temporary directory..."
2480 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2481 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2485 c------------------------------------------------------------------------------
2486 subroutine move_from_tmp(source)
2487 include "DIMENSIONS"
2488 include "COMMON.IOUNITS"
2489 character*(*) source
2492 write (*,*) "Moving ",source(:ilen(source)),
2493 & " from temporary directory to working directory"
2494 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2495 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2498 c------------------------------------------------------------------------------
2499 subroutine random_init(seed)
2501 C Initialize random number generator
2503 implicit real*8 (a-h,o-z)
2504 include 'DIMENSIONS'
2507 logical OKRandom, prng_restart
2509 integer iseed_array(4)
2511 include 'COMMON.IOUNITS'
2512 include 'COMMON.TIME1'
2513 include 'COMMON.THREAD'
2514 include 'COMMON.SBRIDGE'
2515 include 'COMMON.CONTROL'
2516 include 'COMMON.MCM'
2517 include 'COMMON.MAP'
2518 include 'COMMON.HEADER'
2519 include 'COMMON.CSA'
2520 include 'COMMON.CHAIN'
2521 include 'COMMON.MUCA'
2523 include 'COMMON.FFIELD'
2524 include 'COMMON.SETUP'
2525 iseed=-dint(dabs(seed))
2526 if (iseed.eq.0) then
2527 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2528 & 'Random seed undefined. The program will stop.'
2529 write (*,'(/80(1h*)/20x,a/80(1h*))')
2530 & 'Random seed undefined. The program will stop.'
2532 call mpi_finalize(mpi_comm_world,ierr)
2534 stop 'Bad random seed.'
2537 if (fg_rank.eq.0) then
2541 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2542 OKRandom = prng_restart(me,iseed)
2545 tmp=65536.0d0**(4-i)
2546 iseed_array(i) = dint(seed/tmp)
2547 seed=seed-iseed_array(i)*tmp
2550 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2551 & (iseed_array(i),i=1,4)
2552 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2553 & (iseed_array(i),i=1,4)
2554 OKRandom = prng_restart(me,iseed_array)
2557 c r1 = prng_next(me)
2558 r1=ran_number(0.0D0,1.0D0)
2560 & write (iout,*) 'ran_num',r1
2561 if (r1.lt.0.0d0) OKRandom=.false.
2563 if (.not.OKRandom) then
2564 write (iout,*) 'PRNG IS NOT WORKING!!!'
2565 print *,'PRNG IS NOT WORKING!!!'
2568 call mpi_abort(mpi_comm_world,error_msg,ierr)
2571 write (iout,*) 'too many processors for parallel prng'
2572 write (*,*) 'too many processors for parallel prng'
2580 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)