2 implicit real*8 (a-h,o-z)
8 include 'COMMON.CONTROL'
9 include 'COMMON.SBRIDGE'
10 include 'COMMON.IOUNITS'
11 include 'COMMON.SPLITELE'
13 C Read force-field parameters except weights
15 C Read job setup parameters
17 C Read control parameters for energy minimzation if required
18 if (minim) call read_minim
19 C Read MCM control parameters if required
20 if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
21 C Read MD control parameters if reqjuired
22 if (modecalc.eq.12) call read_MDpar
23 C Read MREMD control parameters if required
24 if (modecalc.eq.14) then
28 C Read MUCA control parameters if required
29 if (lmuca) call read_muca
30 C Read CSA control parameters if required (from fort.40 if exists
31 C otherwise from general input file)
32 if (modecalc.eq.8) then
33 inquire (file="fort.40",exist=file_exist)
34 if (.not.file_exist) call csaread
36 cfmc if (modecalc.eq.10) call mcmfread
37 C Read molecule information, molecule geometry, energy-term weights, and
38 C restraints if requested
40 C Print restraint information
42 if (.not. out1file .or. me.eq.king) then
45 &write (iout,'(a,i5,a)') "The following",nhpb-nss,
46 & " distance constraints have been imposed"
48 write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
53 c print *,"Processor",myrank," leaves READRTNS"
56 C-------------------------------------------------------------------------------
57 subroutine read_control
61 implicit real*8 (a-h,o-z)
65 logical OKRandom, prng_restart
68 include 'COMMON.IOUNITS'
69 include 'COMMON.TIME1'
70 include 'COMMON.THREAD'
71 include 'COMMON.SBRIDGE'
72 include 'COMMON.CONTROL'
75 include 'COMMON.HEADER'
77 include 'COMMON.CHAIN'
80 include 'COMMON.FFIELD'
81 include 'COMMON.INTERACT'
82 include 'COMMON.SETUP'
83 include 'COMMON.SPLITELE'
84 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
85 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
87 character*320 controlcard
92 read (INP,'(a)') titel
93 call card_concat(controlcard)
94 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
95 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
96 call reada(controlcard,'SEED',seed,0.0D0)
97 call random_init(seed)
98 C Set up the time limit (caution! The time must be input in minutes!)
99 read_cart=index(controlcard,'READ_CART').gt.0
100 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
101 call readi(controlcard,'SYM',symetr,1)
102 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
103 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
104 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
105 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
106 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
107 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
108 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
109 call reada(controlcard,'DRMS',drms,0.1D0)
110 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
111 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
112 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
113 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
114 write (iout,'(a,f10.1)')'DRMS = ',drms
115 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
116 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
118 call readi(controlcard,'NZ_START',nz_start,0)
119 call readi(controlcard,'NZ_END',nz_end,0)
120 call readi(controlcard,'IZ_SC',iz_sc,0)
122 safety = 60.0d0*safety
125 call reada(controlcard,"T_BATH",t_bath,300.0d0)
126 minim=(index(controlcard,'MINIMIZE').gt.0)
127 dccart=(index(controlcard,'CART').gt.0)
128 overlapsc=(index(controlcard,'OVERLAP').gt.0)
129 overlapsc=.not.overlapsc
130 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
131 searchsc=.not.searchsc
132 sideadd=(index(controlcard,'SIDEADD').gt.0)
133 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
134 outpdb=(index(controlcard,'PDBOUT').gt.0)
135 outmol2=(index(controlcard,'MOL2OUT').gt.0)
136 pdbref=(index(controlcard,'PDBREF').gt.0)
137 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
138 indpdb=index(controlcard,'PDBSTART')
139 extconf=(index(controlcard,'EXTCONF').gt.0)
140 call readi(controlcard,'IPRINT',iprint,0)
141 call readi(controlcard,'MAXGEN',maxgen,10000)
142 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
143 call readi(controlcard,"KDIAG",kdiag,0)
144 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
145 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
146 & write (iout,*) "RESCALE_MODE",rescale_mode
147 split_ene=index(controlcard,'SPLIT_ENE').gt.0
148 if (index(controlcard,'REGULAR').gt.0.0D0) then
149 call reada(controlcard,'WEIDIS',weidis,0.1D0)
153 if (index(controlcard,'CHECKGRAD').gt.0) then
155 if (index(controlcard,'CART').gt.0) then
157 elseif (index(controlcard,'CARINT').gt.0) then
162 elseif (index(controlcard,'THREAD').gt.0) then
164 call readi(controlcard,'THREAD',nthread,0)
165 if (nthread.gt.0) then
166 call reada(controlcard,'WEIDIS',weidis,0.1D0)
169 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
170 stop 'Error termination in Read_Control.'
172 else if (index(controlcard,'MCMA').gt.0) then
174 else if (index(controlcard,'MCEE').gt.0) then
176 else if (index(controlcard,'MULTCONF').gt.0) then
178 else if (index(controlcard,'MAP').gt.0) then
180 call readi(controlcard,'MAP',nmap,0)
181 else if (index(controlcard,'CSA').gt.0) then
183 crc else if (index(controlcard,'ZSCORE').gt.0) then
185 crc ZSCORE is rm from UNRES, modecalc=9 is available
188 cfcm else if (index(controlcard,'MCMF').gt.0) then
190 else if (index(controlcard,'SOFTREG').gt.0) then
192 else if (index(controlcard,'CHECK_BOND').gt.0) then
194 else if (index(controlcard,'TEST').gt.0) then
196 else if (index(controlcard,'MD').gt.0) then
198 else if (index(controlcard,'RE ').gt.0) then
202 lmuca=index(controlcard,'MUCA').gt.0
203 call readi(controlcard,'MUCADYN',mucadyn,0)
204 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
205 if (lmuca .and. (me.eq.king .or. .not.out1file ))
207 write (iout,*) 'MUCADYN=',mucadyn
208 write (iout,*) 'MUCASMOOTH=',muca_smooth
211 iscode=index(controlcard,'ONE_LETTER')
212 indphi=index(controlcard,'PHI')
213 indback=index(controlcard,'BACK')
214 iranconf=index(controlcard,'RAND_CONF')
215 i2ndstr=index(controlcard,'USE_SEC_PRED')
216 gradout=index(controlcard,'GRADOUT').gt.0
217 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
218 C DISTCHAINMAX become obsolete for periodic boundry condition
219 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
220 C Reading the dimensions of box in x,y,z coordinates
221 call reada(controlcard,'BOXX',boxxsize,100.0d0)
222 call reada(controlcard,'BOXY',boxysize,100.0d0)
223 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
224 c Cutoff range for interactions
225 call reada(controlcard,"R_CUT",r_cut,15.0d0)
226 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
227 if (me.eq.king .or. .not.out1file )
228 & write (iout,*) "DISTCHAINMAX",distchainmax
230 if(me.eq.king.or..not.out1file)
231 & write (iout,'(2a)') diagmeth(kdiag),
232 & ' routine used to diagonalize matrices.'
235 c--------------------------------------------------------------------------
236 subroutine read_REMDpar
240 implicit real*8 (a-h,o-z)
242 include 'COMMON.IOUNITS'
243 include 'COMMON.TIME1'
246 include 'COMMON.LANGEVIN'
248 include 'COMMON.LANGEVIN.lang0'
250 include 'COMMON.INTERACT'
251 include 'COMMON.NAMES'
253 include 'COMMON.REMD'
254 include 'COMMON.CONTROL'
255 include 'COMMON.SETUP'
257 character*320 controlcard
258 character*3200 controlcard1
259 integer iremd_m_total
261 if(me.eq.king.or..not.out1file)
262 & write (iout,*) "REMD setup"
264 call card_concat(controlcard)
265 call readi(controlcard,"NREP",nrep,3)
266 call readi(controlcard,"NSTEX",nstex,1000)
267 call reada(controlcard,"RETMIN",retmin,10.0d0)
268 call reada(controlcard,"RETMAX",retmax,1000.0d0)
269 mremdsync=(index(controlcard,'SYNC').gt.0)
270 call readi(controlcard,"NSYN",i_sync_step,100)
271 restart1file=(index(controlcard,'REST1FILE').gt.0)
272 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
273 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
274 if(max_cache_traj_use.gt.max_cache_traj)
275 & max_cache_traj_use=max_cache_traj
276 if(me.eq.king.or..not.out1file) then
277 cd if (traj1file) then
278 crc caching is in testing - NTWX is not ignored
279 cd write (iout,*) "NTWX value is ignored"
280 cd write (iout,*) " trajectory is stored to one file by master"
281 cd write (iout,*) " before exchange at NSTEX intervals"
283 write (iout,*) "NREP= ",nrep
284 write (iout,*) "NSTEX= ",nstex
285 write (iout,*) "SYNC= ",mremdsync
286 write (iout,*) "NSYN= ",i_sync_step
287 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
290 if (index(controlcard,'TLIST').gt.0) then
292 call card_concat(controlcard1)
293 read(controlcard1,*) (remd_t(i),i=1,nrep)
294 if(me.eq.king.or..not.out1file)
295 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
298 if (index(controlcard,'MLIST').gt.0) then
300 call card_concat(controlcard1)
301 read(controlcard1,*) (remd_m(i),i=1,nrep)
302 if(me.eq.king.or..not.out1file) then
303 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
306 iremd_m_total=iremd_m_total+remd_m(i)
308 write (iout,*) 'Total number of replicas ',iremd_m_total
311 if(me.eq.king.or..not.out1file)
312 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
315 c--------------------------------------------------------------------------
316 subroutine read_MDpar
320 implicit real*8 (a-h,o-z)
322 include 'COMMON.IOUNITS'
323 include 'COMMON.TIME1'
326 include 'COMMON.LANGEVIN'
328 include 'COMMON.LANGEVIN.lang0'
330 include 'COMMON.INTERACT'
331 include 'COMMON.NAMES'
333 include 'COMMON.SETUP'
334 include 'COMMON.CONTROL'
335 include 'COMMON.SPLITELE'
337 character*320 controlcard
339 call card_concat(controlcard)
340 call readi(controlcard,"NSTEP",n_timestep,1000000)
341 call readi(controlcard,"NTWE",ntwe,100)
342 call readi(controlcard,"NTWX",ntwx,1000)
343 call reada(controlcard,"DT",d_time,1.0d-1)
344 call reada(controlcard,"DVMAX",dvmax,2.0d1)
345 call reada(controlcard,"DAMAX",damax,1.0d1)
346 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
347 call readi(controlcard,"LANG",lang,0)
348 RESPA = index(controlcard,"RESPA") .gt. 0
349 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
350 ntime_split0=ntime_split
351 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
352 ntime_split0=ntime_split
353 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
354 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
355 rest = index(controlcard,"REST").gt.0
356 tbf = index(controlcard,"TBF").gt.0
357 usampl = index(controlcard,"USAMPL").gt.0
358 mdpdb = index(controlcard,"MDPDB").gt.0
359 call reada(controlcard,"T_BATH",t_bath,300.0d0)
360 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
361 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
362 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
363 if (count_reset_moment.eq.0) count_reset_moment=1000000000
364 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
365 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
366 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
367 if (count_reset_vel.eq.0) count_reset_vel=1000000000
368 large = index(controlcard,"LARGE").gt.0
369 print_compon = index(controlcard,"PRINT_COMPON").gt.0
370 rattle = index(controlcard,"RATTLE").gt.0
371 c if performing umbrella sampling, fragments constrained are read from the fragment file
377 if(me.eq.king.or..not.out1file) then
379 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
381 write (iout,'(a)') "The units are:"
382 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
383 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
384 & " acceleration: angstrom/(48.9 fs)**2"
385 write (iout,'(a)') "energy: kcal/mol, temperature: K"
387 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
388 write (iout,'(a60,f10.5,a)')
389 & "Initial time step of numerical integration:",d_time,
391 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
393 write (iout,'(2a,i4,a)')
394 & "A-MTS algorithm used; initial time step for fast-varying",
395 & " short-range forces split into",ntime_split," steps."
396 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
397 & r_cut," lambda",rlamb
399 write (iout,'(2a,f10.5)')
400 & "Maximum acceleration threshold to reduce the time step",
401 & "/increase split number:",damax
402 write (iout,'(2a,f10.5)')
403 & "Maximum predicted energy drift to reduce the timestep",
404 & "/increase split number:",edriftmax
405 write (iout,'(a60,f10.5)')
406 & "Maximum velocity threshold to reduce velocities:",dvmax
407 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
408 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
409 if (rattle) write (iout,'(a60)')
410 & "Rattle algorithm used to constrain the virtual bonds"
414 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
415 call reada(controlcard,"RWAT",rwat,1.4d0)
416 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
417 surfarea=index(controlcard,"SURFAREA").gt.0
418 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
419 if(me.eq.king.or..not.out1file)then
420 write (iout,'(/a,$)') "Langevin dynamics calculation"
423 & " with direct integration of Langevin equations"
424 else if (lang.eq.2) then
425 write (iout,'(a/)') " with TINKER stochasic MD integrator"
426 else if (lang.eq.3) then
427 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
428 else if (lang.eq.4) then
429 write (iout,'(a/)') " in overdamped mode"
431 write (iout,'(//a,i5)')
432 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
435 write (iout,'(a60,f10.5)') "Temperature:",t_bath
436 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
437 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
438 write (iout,'(a60,f10.5)')
439 & "Scaling factor of the friction forces:",scal_fric
440 if (surfarea) write (iout,'(2a,i10,a)')
441 & "Friction coefficients will be scaled by solvent-accessible",
442 & " surface area every",reset_fricmat," steps."
444 c Calculate friction coefficients and bounds of stochastic forces
445 eta=6*pi*cPoise*etawat
446 if(me.eq.king.or..not.out1file)
447 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
449 gamp=scal_fric*(pstok+rwat)*eta
450 stdfp=dsqrt(2*Rb*t_bath/d_time)
452 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
453 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
455 if(me.eq.king.or..not.out1file)then
456 write (iout,'(/2a/)')
457 & "Radii of site types and friction coefficients and std's of",
458 & " stochastic forces of fully exposed sites"
459 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
461 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
462 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
466 if(me.eq.king.or..not.out1file)then
467 write (iout,'(a)') "Berendsen bath calculation"
468 write (iout,'(a60,f10.5)') "Temperature:",t_bath
469 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
471 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
472 & count_reset_moment," steps"
474 & write (iout,'(a,i10,a)')
475 & "Velocities will be reset at random every",count_reset_vel,
479 if(me.eq.king.or..not.out1file)
480 & write (iout,'(a31)') "Microcanonical mode calculation"
482 if(me.eq.king.or..not.out1file)then
483 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
485 write(iout,*) "MD running with constraints."
486 write(iout,*) "Equilibration time ", eq_time, " mtus."
487 write(iout,*) "Constraining ", nfrag," fragments."
488 write(iout,*) "Length of each fragment, weight and q0:"
490 write (iout,*) "Set of restraints #",iset
492 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
493 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
495 write(iout,*) "constraints between ", npair, "fragments."
496 write(iout,*) "constraint pairs, weights and q0:"
498 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
499 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
501 write(iout,*) "angle constraints within ", nfrag_back,
502 & "backbone fragments."
503 write(iout,*) "fragment, weights:"
505 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
506 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
507 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
510 iset=mod(kolor,nset)+1
513 if(me.eq.king.or..not.out1file)
514 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
517 c------------------------------------------------------------------------------
520 C Read molecular data.
522 implicit real*8 (a-h,o-z)
528 include 'COMMON.IOUNITS'
531 include 'COMMON.INTERACT'
532 include 'COMMON.LOCAL'
533 include 'COMMON.NAMES'
534 include 'COMMON.CHAIN'
535 include 'COMMON.FFIELD'
536 include 'COMMON.SBRIDGE'
537 include 'COMMON.HEADER'
538 include 'COMMON.CONTROL'
539 include 'COMMON.DBASE'
540 include 'COMMON.THREAD'
541 include 'COMMON.CONTACTS'
542 include 'COMMON.TORCNSTR'
543 include 'COMMON.TIME1'
544 include 'COMMON.BOUNDS'
546 include 'COMMON.SETUP'
547 character*4 sequence(maxres)
549 double precision x(maxvar)
550 character*256 pdbfile
551 character*320 weightcard
552 character*80 weightcard_t,ucase
553 dimension itype_pdb(maxres)
554 common /pizda/ itype_pdb
555 logical seq_comp,fail
556 double precision energia(0:n_ene)
562 C Read weights of the subsequent energy terms.
563 call card_concat(weightcard)
564 call reada(weightcard,'WLONG',wlong,1.0D0)
565 call reada(weightcard,'WSC',wsc,wlong)
566 call reada(weightcard,'WSCP',wscp,wlong)
567 call reada(weightcard,'WELEC',welec,1.0D0)
568 call reada(weightcard,'WVDWPP',wvdwpp,welec)
569 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
570 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
571 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
572 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
573 call reada(weightcard,'WTURN3',wturn3,1.0D0)
574 call reada(weightcard,'WTURN4',wturn4,1.0D0)
575 call reada(weightcard,'WTURN6',wturn6,1.0D0)
576 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
577 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
578 call reada(weightcard,'WBOND',wbond,1.0D0)
579 call reada(weightcard,'WTOR',wtor,1.0D0)
580 call reada(weightcard,'WTORD',wtor_d,1.0D0)
581 call reada(weightcard,'WANG',wang,1.0D0)
582 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
583 call reada(weightcard,'SCAL14',scal14,0.4D0)
584 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
585 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
586 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
587 call reada(weightcard,'TEMP0',temp0,300.0d0)
588 if (index(weightcard,'SOFT').gt.0) ipot=6
589 C 12/1/95 Added weight for the multi-body term WCORR
590 call reada(weightcard,'WCORRH',wcorr,1.0D0)
591 if (wcorr4.gt.0.0d0) wcorr=wcorr4
611 if(me.eq.king.or..not.out1file)
612 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
613 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
615 10 format (/'Energy-term weights (unscaled):'//
616 & 'WSCC= ',f10.6,' (SC-SC)'/
617 & 'WSCP= ',f10.6,' (SC-p)'/
618 & 'WELEC= ',f10.6,' (p-p electr)'/
619 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
620 & 'WBOND= ',f10.6,' (stretching)'/
621 & 'WANG= ',f10.6,' (bending)'/
622 & 'WSCLOC= ',f10.6,' (SC local)'/
623 & 'WTOR= ',f10.6,' (torsional)'/
624 & 'WTORD= ',f10.6,' (double torsional)'/
625 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
626 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
627 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
628 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
629 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
630 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
631 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
632 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
633 & 'WTURN6= ',f10.6,' (turns, 6th order)')
634 if(me.eq.king.or..not.out1file)then
635 if (wcorr4.gt.0.0d0) then
636 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
637 & 'between contact pairs of peptide groups'
638 write (iout,'(2(a,f5.3/))')
639 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
640 & 'Range of quenching the correlation terms:',2*delt_corr
641 else if (wcorr.gt.0.0d0) then
642 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
643 & 'between contact pairs of peptide groups'
645 write (iout,'(a,f8.3)')
646 & 'Scaling factor of 1,4 SC-p interactions:',scal14
647 write (iout,'(a,f8.3)')
648 & 'General scaling factor of SC-p interactions:',scalscp
650 r0_corr=cutoff_corr-delt_corr
652 aad(i,1)=scalscp*aad(i,1)
653 aad(i,2)=scalscp*aad(i,2)
654 bad(i,1)=scalscp*bad(i,1)
655 bad(i,2)=scalscp*bad(i,2)
657 call rescale_weights(t_bath)
658 if(me.eq.king.or..not.out1file)
659 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
660 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
662 22 format (/'Energy-term weights (scaled):'//
663 & 'WSCC= ',f10.6,' (SC-SC)'/
664 & 'WSCP= ',f10.6,' (SC-p)'/
665 & 'WELEC= ',f10.6,' (p-p electr)'/
666 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
667 & 'WBOND= ',f10.6,' (stretching)'/
668 & 'WANG= ',f10.6,' (bending)'/
669 & 'WSCLOC= ',f10.6,' (SC local)'/
670 & 'WTOR= ',f10.6,' (torsional)'/
671 & 'WTORD= ',f10.6,' (double torsional)'/
672 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
673 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
674 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
675 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
676 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
677 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
678 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
679 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
680 & 'WTURN6= ',f10.6,' (turns, 6th order)')
681 if(me.eq.king.or..not.out1file)
682 & write (iout,*) "Reference temperature for weights calculation:",
684 call reada(weightcard,"D0CM",d0cm,3.78d0)
685 call reada(weightcard,"AKCM",akcm,15.1d0)
686 call reada(weightcard,"AKTH",akth,11.0d0)
687 call reada(weightcard,"AKCT",akct,12.0d0)
688 call reada(weightcard,"V1SS",v1ss,-1.08d0)
689 call reada(weightcard,"V2SS",v2ss,7.61d0)
690 call reada(weightcard,"V3SS",v3ss,13.7d0)
691 call reada(weightcard,"EBR",ebr,-5.50D0)
692 if(me.eq.king.or..not.out1file) then
693 write (iout,*) "Parameters of the SS-bond potential:"
694 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
696 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
697 write (iout,*) "EBR",ebr
698 c print *,'indpdb=',indpdb,' pdbref=',pdbref
700 if (indpdb.gt.0 .or. pdbref) then
701 read(inp,'(a)') pdbfile
702 if(me.eq.king.or..not.out1file)
703 & write (iout,'(2a)') 'PDB data will be read from file ',
704 & pdbfile(:ilen(pdbfile))
705 open(ipdbin,file=pdbfile,status='old',err=33)
707 33 write (iout,'(a)') 'Error opening PDB file.'
710 c print *,'Begin reading pdb data'
712 c print *,'Finished reading pdb data'
713 if(me.eq.king.or..not.out1file)
714 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
715 & ' nstart_sup=',nstart_sup
717 itype_pdb(i)=itype(i)
721 nct=nstart_sup+nsup-1
722 call contact(.false.,ncont_ref,icont_ref,co)
725 if(me.eq.king.or..not.out1file)
726 & write(iout,*)'Adding sidechains'
730 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
733 do while (fail.and.nsi.le.maxsi)
734 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
737 if(fail) write(iout,*)'Adding sidechain failed for res ',
738 & i,' after ',nsi,' trials'
743 if (indpdb.eq.0) then
744 C Read sequence if not taken from the pdb file.
746 c print *,'nres=',nres
747 if (iscode.gt.0) then
748 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
750 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
752 C Convert sequence to numeric code
754 itype(i)=rescode(i,sequence(i),iscode)
756 C Assign initial virtual bond lengths
762 vbld(i+nres)=dsc(iabs(itype(i)))
763 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
764 c write (iout,*) "i",i," itype",itype(i),
765 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
769 c print '(20i4)',(itype(i),i=1,nres)
772 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
774 if (itype(i).eq.ntyp1) then
778 else if (iabs(itype(i+1)).ne.20) then
780 else if (iabs(itype(i)).ne.20) then
787 if(me.eq.king.or..not.out1file)then
788 write (iout,*) "ITEL"
790 write (iout,*) i,itype(i),itel(i)
792 print *,'Call Read_Bridge.'
795 C 8/13/98 Set limits to generating the dihedral angles
800 read (inp,*) ndih_constr
801 if (ndih_constr.gt.0) then
803 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
804 if(me.eq.king.or..not.out1file)then
806 & 'There are',ndih_constr,' constraints on phi angles.'
808 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
812 phi0(i)=deg2rad*phi0(i)
813 drange(i)=deg2rad*drange(i)
815 if(me.eq.king.or..not.out1file)
816 & write (iout,*) 'FTORS',ftors
819 phibound(1,ii) = phi0(i)-drange(i)
820 phibound(2,ii) = phi0(i)+drange(i)
827 write (iout,'(a)') 'Boundaries in phi angle sampling:'
829 write (iout,'(a3,i5,2f10.1)')
830 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
836 cd print *,'NNT=',NNT,' NCT=',NCT
837 if (itype(1).eq.ntyp1) nnt=2
838 if (itype(nres).eq.ntyp1) nct=nct-1
840 if(me.eq.king.or..not.out1file)
841 & write (iout,'(a,i3)') 'nsup=',nsup
843 if (nsup.le.(nct-nnt+1)) then
844 do i=0,nct-nnt+1-nsup
845 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
851 & 'Error - sequences to be superposed do not match.'
854 do i=0,nsup-(nct-nnt+1)
855 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
857 nstart_sup=nstart_sup+i
863 & 'Error - sequences to be superposed do not match.'
866 if (nsup.eq.0) nsup=nct-nnt
867 if (nstart_sup.eq.0) nstart_sup=nnt
868 if (nstart_seq.eq.0) nstart_seq=nnt
869 if(me.eq.king.or..not.out1file)
870 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
871 & ' nstart_seq=',nstart_seq
873 c--- Zscore rms -------
874 if (nz_start.eq.0) nz_start=nnt
875 if (nz_end.eq.0 .and. nsup.gt.0) then
877 else if (nz_end.eq.0) then
880 if(me.eq.king.or..not.out1file)then
881 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
882 write (iout,*) 'IZ_SC=',iz_sc
884 c----------------------
887 if (.not.pdbref) then
888 call read_angles(inp,*38)
890 38 write (iout,'(a)') 'Error reading reference structure.'
892 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
893 stop 'Error reading reference structure'
897 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
907 call contact(.true.,ncont_ref,icont_ref,co)
909 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
911 if (constr_dist.gt.0) call read_dist_constr
912 write (iout,*) "After read_dist_constr nhpb",nhpb
914 if(me.eq.king.or..not.out1file)
915 & write (iout,*) 'Contact order:',co
917 if(me.eq.king.or..not.out1file)
918 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
921 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
923 if(me.eq.king.or..not.out1file)
924 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
925 & icont_ref(1,i),' ',
926 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
930 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
931 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
932 & modecalc.ne.10) then
933 C If input structure hasn't been supplied from the PDB file read or generate
935 if (iranconf.eq.0 .and. .not. extconf) then
936 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
937 & write (iout,'(a)') 'Initial geometry will be read in.'
939 read(inp,'(8f10.5)',end=36,err=36)
940 & ((c(l,k),l=1,3),k=1,nres),
941 & ((c(l,k+nres),l=1,3),k=nnt,nct)
942 write (iout,*) "Exit READ_CART"
943 write (iout,'(8f10.5)')
944 & ((c(l,k),l=1,3),k=1,nres),
945 & ((c(l,k+nres),l=1,3),k=nnt,nct)
946 call int_from_cart1(.true.)
947 write (iout,*) "Finish INT_TO_CART"
950 dc(j,i)=c(j,i+1)-c(j,i)
951 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
955 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
957 dc(j,i+nres)=c(j,i+nres)-c(j,i)
958 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
964 call read_angles(inp,*36)
967 36 write (iout,'(a)') 'Error reading angle file.'
969 call mpi_finalize( MPI_COMM_WORLD,IERR )
971 stop 'Error reading angle file.'
973 else if (extconf) then
974 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
975 & write (iout,'(a)') 'Extended chain initial geometry.'
977 theta(i)=90d0*deg2rad
983 alph(i)=110d0*deg2rad
986 omeg(i)=-120d0*deg2rad
987 if (itype(i).le.0) omeg(i)=-omeg(i)
990 if(me.eq.king.or..not.out1file)
991 & write (iout,'(a)') 'Random-generated initial geometry.'
995 if (me.eq.king .or. fg_rank.eq.0 .and. (
996 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1000 call gen_rand_conf(itmp,*30)
1002 30 write (iout,*) 'Failed to generate random conformation',
1003 & ', itrial=',itrial
1004 write (*,*) 'Processor:',me,
1005 & ' Failed to generate random conformation',
1015 write (iout,'(a,i3,a)') 'Processor:',me,
1016 & ' error in generating random conformation.'
1017 write (*,'(a,i3,a)') 'Processor:',me,
1018 & ' error in generating random conformation.'
1021 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1027 call gen_rand_conf(itmp,*30)
1029 30 write (iout,*) 'Failed to generate random conformation',
1030 & ', itrial=',itrial
1031 write (*,*) 'Failed to generate random conformation',
1032 & ', itrial=',itrial
1034 write (iout,'(a,i3,a)') 'Processor:',me,
1035 & ' error in generating random conformation.'
1036 write (*,'(a,i3,a)') 'Processor:',me,
1037 & ' error in generating random conformation.'
1042 elseif (modecalc.eq.4) then
1043 read (inp,'(a)') intinname
1044 open (intin,file=intinname,status='old',err=333)
1045 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1046 & write (iout,'(a)') 'intinname',intinname
1047 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1049 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1051 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1053 stop 'Error opening angle file.'
1057 C Generate distance constraints, if the PDB structure is to be regularized.
1058 if (nthread.gt.0) then
1059 call read_threadbase
1062 if (me.eq.king .or. .not. out1file)
1064 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1065 write (iout,'(/a,i3,a)')
1066 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1067 write (iout,'(20i4)') (iss(i),i=1,ns)
1068 write (iout,'(/a/)') 'Pre-formed links are:'
1074 if (me.eq.king.or..not.out1file)
1075 & write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1076 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1081 if (i2ndstr.gt.0) call secstrp2dihc
1082 c call geom_to_var(nvar,x)
1083 c call etotal(energia(0))
1084 c call enerprint(energia(0))
1085 c call briefout(0,etot)
1087 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1088 cd write (iout,'(a)') 'Variable list:'
1089 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1091 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1092 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1093 & 'Processor',myrank,': end reading molecular data.'
1097 c--------------------------------------------------------------------------
1098 logical function seq_comp(itypea,itypeb,length)
1100 integer length,itypea(length),itypeb(length)
1103 if (itypea(i).ne.itypeb(i)) then
1111 c-----------------------------------------------------------------------------
1112 subroutine read_bridge
1113 C Read information about disulfide bridges.
1114 implicit real*8 (a-h,o-z)
1115 include 'DIMENSIONS'
1119 include 'COMMON.IOUNITS'
1120 include 'COMMON.GEO'
1121 include 'COMMON.VAR'
1122 include 'COMMON.INTERACT'
1123 include 'COMMON.LOCAL'
1124 include 'COMMON.NAMES'
1125 include 'COMMON.CHAIN'
1126 include 'COMMON.FFIELD'
1127 include 'COMMON.SBRIDGE'
1128 include 'COMMON.HEADER'
1129 include 'COMMON.CONTROL'
1130 include 'COMMON.DBASE'
1131 include 'COMMON.THREAD'
1132 include 'COMMON.TIME1'
1133 include 'COMMON.SETUP'
1134 C Read bridging residues.
1135 read (inp,*) ns,(iss(i),i=1,ns)
1137 if(me.eq.king.or..not.out1file)
1138 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1139 C Check whether the specified bridging residues are cystines.
1141 if (itype(iss(i)).ne.1) then
1142 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1143 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
1144 & ' can form a disulfide bridge?!!!'
1145 write (*,'(2a,i3,a)')
1146 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
1147 & ' can form a disulfide bridge?!!!'
1149 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1154 C Read preformed bridges.
1156 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1157 write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1160 C Check if the residues involved in bridges are in the specified list of
1161 C bridging residues.
1164 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1165 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1166 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1167 & ' contains residues present in other pairs.'
1168 write (*,'(a,i3,a)') 'Disulfide pair',i,
1169 & ' contains residues present in other pairs.'
1171 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1177 if (ihpb(i).eq.iss(j)) goto 10
1179 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1182 if (jhpb(i).eq.iss(j)) goto 20
1184 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1190 ihpb(i)=ihpb(i)+nres
1191 jhpb(i)=jhpb(i)+nres
1197 c----------------------------------------------------------------------------
1198 subroutine read_x(kanal,*)
1199 implicit real*8 (a-h,o-z)
1200 include 'DIMENSIONS'
1201 include 'COMMON.GEO'
1202 include 'COMMON.VAR'
1203 include 'COMMON.CHAIN'
1204 include 'COMMON.IOUNITS'
1205 include 'COMMON.CONTROL'
1206 include 'COMMON.LOCAL'
1207 include 'COMMON.INTERACT'
1208 c Read coordinates from input
1210 read(kanal,'(8f10.5)',end=10,err=10)
1211 & ((c(l,k),l=1,3),k=1,nres),
1212 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1215 c(j,2*nres)=c(j,nres)
1217 call int_from_cart1(.false.)
1220 dc(j,i)=c(j,i+1)-c(j,i)
1221 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1225 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1227 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1228 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1236 c----------------------------------------------------------------------------
1237 subroutine read_threadbase
1238 implicit real*8 (a-h,o-z)
1239 include 'DIMENSIONS'
1240 include 'COMMON.IOUNITS'
1241 include 'COMMON.GEO'
1242 include 'COMMON.VAR'
1243 include 'COMMON.INTERACT'
1244 include 'COMMON.LOCAL'
1245 include 'COMMON.NAMES'
1246 include 'COMMON.CHAIN'
1247 include 'COMMON.FFIELD'
1248 include 'COMMON.SBRIDGE'
1249 include 'COMMON.HEADER'
1250 include 'COMMON.CONTROL'
1251 include 'COMMON.DBASE'
1252 include 'COMMON.THREAD'
1253 include 'COMMON.TIME1'
1254 C Read pattern database for threading.
1255 read (icbase,*) nseq
1257 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1258 & nres_base(2,i),nres_base(3,i)
1259 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1261 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1262 c & nres_base(2,i),nres_base(3,i)
1263 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1267 if (weidis.eq.0.0D0) weidis=0.1D0
1276 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1277 write (iout,'(a,i5)') 'nexcl: ',nexcl
1278 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1281 c------------------------------------------------------------------------------
1282 subroutine setup_var
1283 implicit real*8 (a-h,o-z)
1284 include 'DIMENSIONS'
1285 include 'COMMON.IOUNITS'
1286 include 'COMMON.GEO'
1287 include 'COMMON.VAR'
1288 include 'COMMON.INTERACT'
1289 include 'COMMON.LOCAL'
1290 include 'COMMON.NAMES'
1291 include 'COMMON.CHAIN'
1292 include 'COMMON.FFIELD'
1293 include 'COMMON.SBRIDGE'
1294 include 'COMMON.HEADER'
1295 include 'COMMON.CONTROL'
1296 include 'COMMON.DBASE'
1297 include 'COMMON.THREAD'
1298 include 'COMMON.TIME1'
1299 C Set up variable list.
1305 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1307 ialph(i,1)=nvar+nside
1311 if (indphi.gt.0) then
1313 else if (indback.gt.0) then
1318 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1321 c----------------------------------------------------------------------------
1322 subroutine gen_dist_constr
1323 C Generate CA distance constraints.
1324 implicit real*8 (a-h,o-z)
1325 include 'DIMENSIONS'
1326 include 'COMMON.IOUNITS'
1327 include 'COMMON.GEO'
1328 include 'COMMON.VAR'
1329 include 'COMMON.INTERACT'
1330 include 'COMMON.LOCAL'
1331 include 'COMMON.NAMES'
1332 include 'COMMON.CHAIN'
1333 include 'COMMON.FFIELD'
1334 include 'COMMON.SBRIDGE'
1335 include 'COMMON.HEADER'
1336 include 'COMMON.CONTROL'
1337 include 'COMMON.DBASE'
1338 include 'COMMON.THREAD'
1339 include 'COMMON.TIME1'
1340 dimension itype_pdb(maxres)
1341 common /pizda/ itype_pdb
1343 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1344 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1345 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1347 do i=nstart_sup,nstart_sup+nsup-1
1348 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1349 cd & ' seq_pdb', restyp(itype_pdb(i))
1350 do j=i+2,nstart_sup+nsup-1
1352 ihpb(nhpb)=i+nstart_seq-nstart_sup
1353 jhpb(nhpb)=j+nstart_seq-nstart_sup
1355 dhpb(nhpb)=dist(i,j)
1358 cd write (iout,'(a)') 'Distance constraints:'
1363 cd if (ii.gt.nres) then
1368 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1369 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1370 cd & dhpb(i),forcon(i)
1374 c----------------------------------------------------------------------------
1376 implicit real*8 (a-h,o-z)
1377 include 'DIMENSIONS'
1378 include 'COMMON.MAP'
1379 include 'COMMON.IOUNITS'
1380 character*3 angid(4) /'THE','PHI','ALP','OME'/
1381 character*80 mapcard,ucase
1383 read (inp,'(a)') mapcard
1384 mapcard=ucase(mapcard)
1385 if (index(mapcard,'PHI').gt.0) then
1387 else if (index(mapcard,'THE').gt.0) then
1389 else if (index(mapcard,'ALP').gt.0) then
1391 else if (index(mapcard,'OME').gt.0) then
1394 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1395 stop 'Error - illegal variable spec in MAP card.'
1397 call readi (mapcard,'RES1',res1(imap),0)
1398 call readi (mapcard,'RES2',res2(imap),0)
1399 if (res1(imap).eq.0) then
1400 res1(imap)=res2(imap)
1401 else if (res2(imap).eq.0) then
1402 res2(imap)=res1(imap)
1404 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1406 & 'Error - illegal definition of variable group in MAP.'
1407 stop 'Error - illegal definition of variable group in MAP.'
1409 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1410 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1411 call readi(mapcard,'NSTEP',nstep(imap),0)
1412 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1414 & 'Illegal boundary and/or step size specification in MAP.'
1415 stop 'Illegal boundary and/or step size specification in MAP.'
1420 c----------------------------------------------------------------------------
1422 implicit real*8 (a-h,o-z)
1423 include 'DIMENSIONS'
1424 include 'COMMON.IOUNITS'
1425 include 'COMMON.GEO'
1426 include 'COMMON.CSA'
1427 include 'COMMON.BANK'
1428 include 'COMMON.CONTROL'
1430 character*620 mcmcard
1431 call card_concat(mcmcard)
1433 call readi(mcmcard,'NCONF',nconf,50)
1434 call readi(mcmcard,'NADD',nadd,0)
1435 call readi(mcmcard,'JSTART',jstart,1)
1436 call readi(mcmcard,'JEND',jend,1)
1437 call readi(mcmcard,'NSTMAX',nstmax,500000)
1438 call readi(mcmcard,'N0',n0,1)
1439 call readi(mcmcard,'N1',n1,6)
1440 call readi(mcmcard,'N2',n2,4)
1441 call readi(mcmcard,'N3',n3,0)
1442 call readi(mcmcard,'N4',n4,0)
1443 call readi(mcmcard,'N5',n5,0)
1444 call readi(mcmcard,'N6',n6,10)
1445 call readi(mcmcard,'N7',n7,0)
1446 call readi(mcmcard,'N8',n8,0)
1447 call readi(mcmcard,'N9',n9,0)
1448 call readi(mcmcard,'N14',n14,0)
1449 call readi(mcmcard,'N15',n15,0)
1450 call readi(mcmcard,'N16',n16,0)
1451 call readi(mcmcard,'N17',n17,0)
1452 call readi(mcmcard,'N18',n18,0)
1454 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1456 call readi(mcmcard,'NDIFF',ndiff,2)
1457 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1458 call readi(mcmcard,'IS1',is1,1)
1459 call readi(mcmcard,'IS2',is2,8)
1460 call readi(mcmcard,'NRAN0',nran0,4)
1461 call readi(mcmcard,'NRAN1',nran1,2)
1462 call readi(mcmcard,'IRR',irr,1)
1463 call readi(mcmcard,'NSEED',nseed,20)
1464 call readi(mcmcard,'NTOTAL',ntotal,10000)
1465 call reada(mcmcard,'CUT1',cut1,2.0d0)
1466 call reada(mcmcard,'CUT2',cut2,5.0d0)
1467 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1468 call readi(mcmcard,'ICMAX',icmax,3)
1469 call readi(mcmcard,'IRESTART',irestart,0)
1470 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1473 call reada(mcmcard,'DELE',dele,20.0d0)
1474 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1475 call readi(mcmcard,'IREF',iref,0)
1476 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1477 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1478 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1479 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1480 write (iout,*) "NCONF_IN",nconf_in
1483 c----------------------------------------------------------------------------
1484 cfmc subroutine mcmfread
1485 cfmc implicit real*8 (a-h,o-z)
1486 cfmc include 'DIMENSIONS'
1487 cfmc include 'COMMON.MCMF'
1488 cfmc include 'COMMON.IOUNITS'
1489 cfmc include 'COMMON.GEO'
1490 cfmc character*80 ucase
1491 cfmc character*620 mcmcard
1492 cfmc call card_concat(mcmcard)
1494 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1495 cfmc write(iout,*)'MAXRANT=',maxrant
1496 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1497 cfmc write(iout,*)'MAXFAM=',maxfam
1498 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1499 cfmc write(iout,*)'NNET1=',nnet1
1500 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1501 cfmc write(iout,*)'NNET2=',nnet2
1502 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1503 cfmc write(iout,*)'NNET3=',nnet3
1504 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1505 cfmc write(iout,*)'ILASTT=',ilastt
1506 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1507 cfmc write(iout,*)'MAXSTR=',maxstr
1508 cfmc maxstr_f=maxstr/maxfam
1509 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1510 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1511 cfmc write(iout,*)'NMCMF=',nmcmf
1512 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1513 cfmc write(iout,*)'IFOCUS=',ifocus
1514 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1515 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1516 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1517 cfmc write(iout,*)'INTPRT=',intprt
1518 cfmc call readi(mcmcard,'IPRT',iprt,100)
1519 cfmc write(iout,*)'IPRT=',iprt
1520 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1521 cfmc write(iout,*)'IMAXTR=',imaxtr
1522 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1523 cfmc write(iout,*)'MAXEVEN=',maxeven
1524 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1525 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1526 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1527 cfmc write(iout,*)'INIMIN=',inimin
1528 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1529 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1530 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1531 cfmc write(iout,*)'NTHREAD=',nthread
1532 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1533 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1534 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1535 cfmc write(iout,*)'MAXPERT=',maxpert
1536 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1537 cfmc write(iout,*)'IRMSD=',irmsd
1538 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1539 cfmc write(iout,*)'DENEMIN=',denemin
1540 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1541 cfmc write(iout,*)'RCUT1S=',rcut1s
1542 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1543 cfmc write(iout,*)'RCUT1E=',rcut1e
1544 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1545 cfmc write(iout,*)'RCUT2S=',rcut2s
1546 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1547 cfmc write(iout,*)'RCUT2E=',rcut2e
1548 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1549 cfmc write(iout,*)'DPERT1=',d_pert1
1550 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1551 cfmc write(iout,*)'DPERT1A=',d_pert1a
1552 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1553 cfmc write(iout,*)'DPERT2=',d_pert2
1554 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1555 cfmc write(iout,*)'DPERT2A=',d_pert2a
1556 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1557 cfmc write(iout,*)'DPERT2B=',d_pert2b
1558 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1559 cfmc write(iout,*)'DPERT2C=',d_pert2c
1560 cfmc d_pert1=deg2rad*d_pert1
1561 cfmc d_pert1a=deg2rad*d_pert1a
1562 cfmc d_pert2=deg2rad*d_pert2
1563 cfmc d_pert2a=deg2rad*d_pert2a
1564 cfmc d_pert2b=deg2rad*d_pert2b
1565 cfmc d_pert2c=deg2rad*d_pert2c
1566 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1567 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1568 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1569 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1570 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1571 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1572 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1573 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1574 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1575 cfmc write(iout,*)'RCUTINI=',rcutini
1576 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1577 cfmc write(iout,*)'GRAT=',grat
1578 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1579 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1583 c----------------------------------------------------------------------------
1585 implicit real*8 (a-h,o-z)
1586 include 'DIMENSIONS'
1587 include 'COMMON.MCM'
1588 include 'COMMON.MCE'
1589 include 'COMMON.IOUNITS'
1591 character*320 mcmcard
1592 call card_concat(mcmcard)
1593 call readi(mcmcard,'MAXACC',maxacc,100)
1594 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1595 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1596 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1597 call readi(mcmcard,'MAXREPM',maxrepm,200)
1598 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1599 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1600 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1601 call reada(mcmcard,'E_UP',e_up,5.0D0)
1602 call reada(mcmcard,'DELTE',delte,0.1D0)
1603 call readi(mcmcard,'NSWEEP',nsweep,5)
1604 call readi(mcmcard,'NSTEPH',nsteph,0)
1605 call readi(mcmcard,'NSTEPC',nstepc,0)
1606 call reada(mcmcard,'TMIN',tmin,298.0D0)
1607 call reada(mcmcard,'TMAX',tmax,298.0D0)
1608 call readi(mcmcard,'NWINDOW',nwindow,0)
1609 call readi(mcmcard,'PRINT_MC',print_mc,0)
1610 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1611 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1612 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1613 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1614 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1615 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1616 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1617 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1618 if (nwindow.gt.0) then
1619 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1621 winlen(i)=winend(i)-winstart(i)+1
1624 if (tmax.lt.tmin) tmax=tmin
1625 if (tmax.eq.tmin) then
1629 if (nstepc.gt.0 .and. nsteph.gt.0) then
1630 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1631 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1633 C Probabilities of different move types
1634 sumpro_type(0)=0.0D0
1635 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1636 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1637 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1638 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1639 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1640 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1641 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1643 print *,'i',i,' sumprotype',sumpro_type(i)
1644 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1645 print *,'i',i,' sumprotype',sumpro_type(i)
1649 c----------------------------------------------------------------------------
1650 subroutine read_minim
1651 implicit real*8 (a-h,o-z)
1652 include 'DIMENSIONS'
1653 include 'COMMON.MINIM'
1654 include 'COMMON.IOUNITS'
1656 character*320 minimcard
1657 call card_concat(minimcard)
1658 call readi(minimcard,'MAXMIN',maxmin,2000)
1659 call readi(minimcard,'MAXFUN',maxfun,5000)
1660 call readi(minimcard,'MINMIN',minmin,maxmin)
1661 call readi(minimcard,'MINFUN',minfun,maxmin)
1662 call reada(minimcard,'TOLF',tolf,1.0D-2)
1663 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1664 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1665 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1666 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1667 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1668 & 'Options in energy minimization:'
1669 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1670 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1671 & 'MinMin:',MinMin,' MinFun:',MinFun,
1672 & ' TolF:',TolF,' RTolF:',RTolF
1675 c----------------------------------------------------------------------------
1676 subroutine read_angles(kanal,*)
1677 implicit real*8 (a-h,o-z)
1678 include 'DIMENSIONS'
1679 include 'COMMON.GEO'
1680 include 'COMMON.VAR'
1681 include 'COMMON.CHAIN'
1682 include 'COMMON.IOUNITS'
1683 include 'COMMON.CONTROL'
1684 c Read angles from input
1686 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1687 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1688 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1689 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1692 c 9/7/01 avoid 180 deg valence angle
1693 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1695 theta(i)=deg2rad*theta(i)
1696 phi(i)=deg2rad*phi(i)
1697 alph(i)=deg2rad*alph(i)
1698 omeg(i)=deg2rad*omeg(i)
1703 c----------------------------------------------------------------------------
1704 subroutine reada(rekord,lancuch,wartosc,default)
1706 character*(*) rekord,lancuch
1707 double precision wartosc,default
1710 iread=index(rekord,lancuch)
1711 if (iread.eq.0) then
1715 iread=iread+ilen(lancuch)+1
1716 read (rekord(iread:),*,err=10,end=10) wartosc
1721 c----------------------------------------------------------------------------
1722 subroutine readi(rekord,lancuch,wartosc,default)
1724 character*(*) rekord,lancuch
1725 integer wartosc,default
1728 iread=index(rekord,lancuch)
1729 if (iread.eq.0) then
1733 iread=iread+ilen(lancuch)+1
1734 read (rekord(iread:),*,err=10,end=10) wartosc
1739 c----------------------------------------------------------------------------
1740 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1743 integer tablica(dim),default
1744 character*(*) rekord,lancuch
1751 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1752 if (iread.eq.0) return
1753 iread=iread+ilen(lancuch)+1
1754 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1757 c----------------------------------------------------------------------------
1758 subroutine multreada(rekord,lancuch,tablica,dim,default)
1761 double precision tablica(dim),default
1762 character*(*) rekord,lancuch
1769 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1770 if (iread.eq.0) return
1771 iread=iread+ilen(lancuch)+1
1772 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1775 c----------------------------------------------------------------------------
1776 subroutine openunits
1777 implicit real*8 (a-h,o-z)
1778 include 'DIMENSIONS'
1781 character*16 form,nodename
1784 include 'COMMON.SETUP'
1785 include 'COMMON.IOUNITS'
1787 include 'COMMON.CONTROL'
1788 integer lenpre,lenpot,ilen,lentmp
1790 character*3 out1file_text,ucase
1793 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1794 call getenv_loc("PREFIX",prefix)
1796 call getenv_loc("POT",pot)
1797 call getenv_loc("DIRTMP",tmpdir)
1798 call getenv_loc("CURDIR",curdir)
1799 call getenv_loc("OUT1FILE",out1file_text)
1800 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1801 out1file_text=ucase(out1file_text)
1802 if (out1file_text(1:1).eq."Y") then
1805 out1file=fg_rank.gt.0
1810 if (lentmp.gt.0) then
1811 write (*,'(80(1h!))')
1812 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1813 write (*,'(80(1h!))')
1814 write (*,*)"All output files will be on node /tmp directory."
1816 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1817 if (me.eq.king) then
1818 write (*,*) "The master node is ",nodename
1819 else if (fg_rank.eq.0) then
1820 write (*,*) "I am the CG slave node ",nodename
1822 write (*,*) "I am the FG slave node ",nodename
1825 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1826 lenpre = lentmp+lenpre+1
1828 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1829 C Get the names and open the input files
1830 #if defined(WINIFL) || defined(WINPGI)
1831 open(1,file=pref_orig(:ilen(pref_orig))//
1832 & '.inp',status='old',readonly,shared)
1833 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1834 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1835 C Get parameter filenames and open the parameter files.
1836 call getenv_loc('BONDPAR',bondname)
1837 open (ibond,file=bondname,status='old',readonly,shared)
1838 call getenv_loc('THETPAR',thetname)
1839 open (ithep,file=thetname,status='old',readonly,shared)
1840 call getenv_loc('ROTPAR',rotname)
1841 open (irotam,file=rotname,status='old',readonly,shared)
1842 call getenv_loc('TORPAR',torname)
1843 open (itorp,file=torname,status='old',readonly,shared)
1844 call getenv_loc('TORDPAR',tordname)
1845 open (itordp,file=tordname,status='old',readonly,shared)
1846 call getenv_loc('FOURIER',fouriername)
1847 open (ifourier,file=fouriername,status='old',readonly,shared)
1848 call getenv_loc('ELEPAR',elename)
1849 open (ielep,file=elename,status='old',readonly,shared)
1850 call getenv_loc('SIDEPAR',sidename)
1851 open (isidep,file=sidename,status='old',readonly,shared)
1852 #elif (defined CRAY) || (defined AIX)
1853 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1855 c print *,"Processor",myrank," opened file 1"
1856 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1857 c print *,"Processor",myrank," opened file 9"
1858 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1859 C Get parameter filenames and open the parameter files.
1860 call getenv_loc('BONDPAR',bondname)
1861 open (ibond,file=bondname,status='old',action='read')
1862 c print *,"Processor",myrank," opened file IBOND"
1863 call getenv_loc('THETPAR',thetname)
1864 open (ithep,file=thetname,status='old',action='read')
1865 c print *,"Processor",myrank," opened file ITHEP"
1866 call getenv_loc('ROTPAR',rotname)
1867 open (irotam,file=rotname,status='old',action='read')
1868 c print *,"Processor",myrank," opened file IROTAM"
1869 call getenv_loc('TORPAR',torname)
1870 open (itorp,file=torname,status='old',action='read')
1871 c print *,"Processor",myrank," opened file ITORP"
1872 call getenv_loc('TORDPAR',tordname)
1873 open (itordp,file=tordname,status='old',action='read')
1874 c print *,"Processor",myrank," opened file ITORDP"
1875 call getenv_loc('SCCORPAR',sccorname)
1876 open (isccor,file=sccorname,status='old',action='read')
1877 c print *,"Processor",myrank," opened file ISCCOR"
1878 call getenv_loc('FOURIER',fouriername)
1879 open (ifourier,file=fouriername,status='old',action='read')
1880 c print *,"Processor",myrank," opened file IFOURIER"
1881 call getenv_loc('ELEPAR',elename)
1882 open (ielep,file=elename,status='old',action='read')
1883 c print *,"Processor",myrank," opened file IELEP"
1884 call getenv_loc('SIDEPAR',sidename)
1885 open (isidep,file=sidename,status='old',action='read')
1886 c print *,"Processor",myrank," opened file ISIDEP"
1887 c print *,"Processor",myrank," opened parameter files"
1889 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1890 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1891 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1892 C Get parameter filenames and open the parameter files.
1893 call getenv_loc('BONDPAR',bondname)
1894 open (ibond,file=bondname,status='old')
1895 call getenv_loc('THETPAR',thetname)
1896 open (ithep,file=thetname,status='old')
1897 call getenv_loc('ROTPAR',rotname)
1898 open (irotam,file=rotname,status='old')
1899 call getenv_loc('TORPAR',torname)
1900 open (itorp,file=torname,status='old')
1901 call getenv_loc('TORDPAR',tordname)
1902 open (itordp,file=tordname,status='old')
1903 call getenv_loc('SCCORPAR',sccorname)
1904 open (isccor,file=sccorname,status='old')
1905 call getenv_loc('FOURIER',fouriername)
1906 open (ifourier,file=fouriername,status='old')
1907 call getenv_loc('ELEPAR',elename)
1908 open (ielep,file=elename,status='old')
1909 call getenv_loc('SIDEPAR',sidename)
1910 open (isidep,file=sidename,status='old')
1912 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1914 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1915 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1916 C Get parameter filenames and open the parameter files.
1917 call getenv_loc('BONDPAR',bondname)
1918 open (ibond,file=bondname,status='old',readonly)
1919 call getenv_loc('THETPAR',thetname)
1920 open (ithep,file=thetname,status='old',readonly)
1921 call getenv_loc('ROTPAR',rotname)
1922 open (irotam,file=rotname,status='old',readonly)
1923 call getenv_loc('TORPAR',torname)
1924 open (itorp,file=torname,status='old',readonly)
1925 call getenv_loc('TORDPAR',tordname)
1926 open (itordp,file=tordname,status='old',readonly)
1927 call getenv_loc('SCCORPAR',sccorname)
1928 open (isccor,file=sccorname,status='old',readonly)
1930 call getenv_loc('THETPARPDB',thetname_pdb)
1931 print *,"thetname_pdb ",thetname_pdb
1932 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1933 print *,ithep_pdb," opened"
1935 call getenv_loc('FOURIER',fouriername)
1936 open (ifourier,file=fouriername,status='old',readonly)
1937 call getenv_loc('ELEPAR',elename)
1938 open (ielep,file=elename,status='old',readonly)
1939 call getenv_loc('SIDEPAR',sidename)
1940 open (isidep,file=sidename,status='old',readonly)
1942 call getenv_loc('ROTPARPDB',rotname_pdb)
1943 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
1948 C 8/9/01 In the newest version SCp interaction constants are read from a file
1949 C Use -DOLDSCP to use hard-coded constants instead.
1951 call getenv_loc('SCPPAR',scpname)
1952 #if defined(WINIFL) || defined(WINPGI)
1953 open (iscpp,file=scpname,status='old',readonly,shared)
1954 #elif (defined CRAY) || (defined AIX)
1955 open (iscpp,file=scpname,status='old',action='read')
1957 open (iscpp,file=scpname,status='old')
1959 open (iscpp,file=scpname,status='old',readonly)
1962 call getenv_loc('PATTERN',patname)
1963 #if defined(WINIFL) || defined(WINPGI)
1964 open (icbase,file=patname,status='old',readonly,shared)
1965 #elif (defined CRAY) || (defined AIX)
1966 open (icbase,file=patname,status='old',action='read')
1968 open (icbase,file=patname,status='old')
1970 open (icbase,file=patname,status='old',readonly)
1973 C Open output file only for CG processes
1974 c print *,"Processor",myrank," fg_rank",fg_rank
1975 if (fg_rank.eq.0) then
1977 if (nodes.eq.1) then
1980 npos = dlog10(dfloat(nodes-1))+1
1982 if (npos.lt.3) npos=3
1983 write (liczba,'(i1)') npos
1984 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
1986 write (liczba,form) me
1987 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
1988 & liczba(:ilen(liczba))
1989 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1991 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1993 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
1994 & liczba(:ilen(liczba))//'.mol2'
1995 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
1996 & liczba(:ilen(liczba))//'.stat'
1998 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
1999 & //liczba(:ilen(liczba))//'.stat')
2000 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2003 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2004 & liczba(:ilen(liczba))//'.const'
2009 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2010 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2011 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2012 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2013 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2015 & call copy_to_tmp(pref_orig(:ile(pref_orig))//'_'//pot(:lenpot)//
2017 rest2name=prefix(:ilen(prefix))//'.rst'
2019 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2022 #if defined(AIX) || defined(PGI)
2023 if (me.eq.king .or. .not. out1file)
2024 & open(iout,file=outname,status='unknown')
2026 if (fg_rank.gt.0) then
2027 write (liczba,'(i3.3)') myrank/nfgtasks
2028 write (ll,'(bz,i3.3)') fg_rank
2029 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2034 open(igeom,file=intname,status='unknown',position='append')
2035 open(ipdb,file=pdbname,status='unknown')
2036 open(imol2,file=mol2name,status='unknown')
2037 open(istat,file=statname,status='unknown',position='append')
2039 c1out open(iout,file=outname,status='unknown')
2042 if (me.eq.king .or. .not.out1file)
2043 & open(iout,file=outname,status='unknown')
2045 if (fg_rank.gt.0) then
2046 write (liczba,'(i3.3)') myrank/nfgtasks
2047 write (ll,'(bz,i3.3)') fg_rank
2048 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2053 open(igeom,file=intname,status='unknown',access='append')
2054 open(ipdb,file=pdbname,status='unknown')
2055 open(imol2,file=mol2name,status='unknown')
2056 open(istat,file=statname,status='unknown',access='append')
2058 c1out open(iout,file=outname,status='unknown')
2061 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2062 csa_seed=prefix(:lenpre)//'.CSA.seed'
2063 csa_history=prefix(:lenpre)//'.CSA.history'
2064 csa_bank=prefix(:lenpre)//'.CSA.bank'
2065 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2066 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2067 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2068 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2069 csa_int=prefix(:lenpre)//'.int'
2070 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2071 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2072 csa_in=prefix(:lenpre)//'.CSA.in'
2073 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2076 write (iout,'(80(1h-))')
2077 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2078 write (iout,'(80(1h-))')
2079 write (iout,*) "Input file : ",
2080 & pref_orig(:ilen(pref_orig))//'.inp'
2081 write (iout,*) "Output file : ",
2082 & outname(:ilen(outname))
2084 write (iout,*) "Sidechain potential file : ",
2085 & sidename(:ilen(sidename))
2087 write (iout,*) "SCp potential file : ",
2088 & scpname(:ilen(scpname))
2090 write (iout,*) "Electrostatic potential file : ",
2091 & elename(:ilen(elename))
2092 write (iout,*) "Cumulant coefficient file : ",
2093 & fouriername(:ilen(fouriername))
2094 write (iout,*) "Torsional parameter file : ",
2095 & torname(:ilen(torname))
2096 write (iout,*) "Double torsional parameter file : ",
2097 & tordname(:ilen(tordname))
2098 write (iout,*) "SCCOR parameter file : ",
2099 & sccorname(:ilen(sccorname))
2100 write (iout,*) "Bond & inertia constant file : ",
2101 & bondname(:ilen(bondname))
2102 write (iout,*) "Bending parameter file : ",
2103 & thetname(:ilen(thetname))
2104 write (iout,*) "Rotamer parameter file : ",
2105 & rotname(:ilen(rotname))
2106 write (iout,*) "Thetpdb parameter file : ",
2107 & thetname_pdb(:ilen(thetname_pdb))
2108 write (iout,*) "Threading database : ",
2109 & patname(:ilen(patname))
2111 &write (iout,*)" DIRTMP : ",
2113 write (iout,'(80(1h-))')
2117 c----------------------------------------------------------------------------
2118 subroutine card_concat(card)
2119 implicit real*8 (a-h,o-z)
2120 include 'DIMENSIONS'
2121 include 'COMMON.IOUNITS'
2123 character*80 karta,ucase
2125 read (inp,'(a)') karta
2128 do while (karta(80:80).eq.'&')
2129 card=card(:ilen(card)+1)//karta(:79)
2130 read (inp,'(a)') karta
2133 card=card(:ilen(card)+1)//karta
2136 c----------------------------------------------------------------------------------
2138 implicit real*8 (a-h,o-z)
2139 include 'DIMENSIONS'
2140 include 'COMMON.CHAIN'
2141 include 'COMMON.IOUNITS'
2143 open(irest2,file=rest2name,status='unknown')
2144 read(irest2,*) totT,EK,potE,totE,t_bath
2146 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2149 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2152 read (irest2,*) iset
2157 c---------------------------------------------------------------------------------
2158 subroutine read_fragments
2159 implicit real*8 (a-h,o-z)
2160 include 'DIMENSIONS'
2164 include 'COMMON.SETUP'
2165 include 'COMMON.CHAIN'
2166 include 'COMMON.IOUNITS'
2168 include 'COMMON.CONTROL'
2169 read(inp,*) nset,nfrag,npair,nfrag_back
2170 if(me.eq.king.or..not.out1file)
2171 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2172 & " nfrag_back",nfrag_back
2174 read(inp,*) mset(iset)
2176 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2178 if(me.eq.king.or..not.out1file)
2179 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2180 & ifrag(2,i,iset), qinfrag(i,iset)
2183 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2185 if(me.eq.king.or..not.out1file)
2186 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2187 & ipair(2,i,iset), qinpair(i,iset)
2190 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2191 & wfrag_back(3,i,iset),
2192 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2193 if(me.eq.king.or..not.out1file)
2194 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2195 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2200 c-------------------------------------------------------------------------------
2201 subroutine read_dist_constr
2202 implicit real*8 (a-h,o-z)
2203 include 'DIMENSIONS'
2207 include 'COMMON.SETUP'
2208 include 'COMMON.CONTROL'
2209 include 'COMMON.CHAIN'
2210 include 'COMMON.IOUNITS'
2211 include 'COMMON.SBRIDGE'
2212 integer ifrag_(2,100),ipair_(2,100)
2213 double precision wfrag_(100),wpair_(100)
2214 character*500 controlcard
2215 c write (iout,*) "Calling read_dist_constr"
2216 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2218 call card_concat(controlcard)
2219 call readi(controlcard,"NFRAG",nfrag_,0)
2220 call readi(controlcard,"NPAIR",npair_,0)
2221 call readi(controlcard,"NDIST",ndist_,0)
2222 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2223 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2224 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2225 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2226 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2227 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2228 c write (iout,*) "IFRAG"
2230 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2232 c write (iout,*) "IPAIR"
2234 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2238 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2239 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2240 & ifrag_(2,i)=nstart_sup+nsup-1
2241 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2243 if (wfrag_(i).gt.0.0d0) then
2244 do j=ifrag_(1,i),ifrag_(2,i)-1
2245 do k=j+1,ifrag_(2,i)
2246 c write (iout,*) "j",j," k",k
2248 if (constr_dist.eq.1) then
2253 forcon(nhpb)=wfrag_(i)
2254 else if (constr_dist.eq.2) then
2255 if (ddjk.le.dist_cut) then
2260 forcon(nhpb)=wfrag_(i)
2267 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2270 if (.not.out1file .or. me.eq.king)
2271 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2272 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2274 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2275 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2282 if (wpair_(i).gt.0.0d0) then
2290 do j=ifrag_(1,ii),ifrag_(2,ii)
2291 do k=ifrag_(1,jj),ifrag_(2,jj)
2295 forcon(nhpb)=wpair_(i)
2296 dhpb(nhpb)=dist(j,k)
2298 if (.not.out1file .or. me.eq.king)
2299 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2300 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2302 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2303 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2310 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2311 if (forcon(nhpb+1).gt.0.0d0) then
2313 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2315 if (.not.out1file .or. me.eq.king)
2316 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2317 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2319 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2320 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2327 c-------------------------------------------------------------------------------
2329 subroutine flush(iu)
2334 subroutine flush(iu)
2339 c------------------------------------------------------------------------------
2340 subroutine copy_to_tmp(source)
2341 include "DIMENSIONS"
2342 include "COMMON.IOUNITS"
2343 character*(*) source
2344 character* 256 tmpfile
2348 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2349 inquire(file=tmpfile,exist=ex)
2351 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2352 & " to temporary directory..."
2353 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2354 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2358 c------------------------------------------------------------------------------
2359 subroutine move_from_tmp(source)
2360 include "DIMENSIONS"
2361 include "COMMON.IOUNITS"
2362 character*(*) source
2365 write (*,*) "Moving ",source(:ilen(source)),
2366 & " from temporary directory to working directory"
2367 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2368 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2371 c------------------------------------------------------------------------------
2372 subroutine random_init(seed)
2374 C Initialize random number generator
2376 implicit real*8 (a-h,o-z)
2377 include 'DIMENSIONS'
2380 logical OKRandom, prng_restart
2382 integer iseed_array(4)
2384 include 'COMMON.IOUNITS'
2385 include 'COMMON.TIME1'
2386 include 'COMMON.THREAD'
2387 include 'COMMON.SBRIDGE'
2388 include 'COMMON.CONTROL'
2389 include 'COMMON.MCM'
2390 include 'COMMON.MAP'
2391 include 'COMMON.HEADER'
2392 include 'COMMON.CSA'
2393 include 'COMMON.CHAIN'
2394 include 'COMMON.MUCA'
2396 include 'COMMON.FFIELD'
2397 include 'COMMON.SETUP'
2398 iseed=-dint(dabs(seed))
2399 if (iseed.eq.0) then
2400 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2401 & 'Random seed undefined. The program will stop.'
2402 write (*,'(/80(1h*)/20x,a/80(1h*))')
2403 & 'Random seed undefined. The program will stop.'
2405 call mpi_finalize(mpi_comm_world,ierr)
2407 stop 'Bad random seed.'
2410 if (fg_rank.eq.0) then
2414 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2415 OKRandom = prng_restart(me,iseed)
2418 tmp=65536.0d0**(4-i)
2419 iseed_array(i) = dint(seed/tmp)
2420 seed=seed-iseed_array(i)*tmp
2423 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2424 & (iseed_array(i),i=1,4)
2425 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2426 & (iseed_array(i),i=1,4)
2427 OKRandom = prng_restart(me,iseed_array)
2430 c r1 = prng_next(me)
2431 r1=ran_number(0.0D0,1.0D0)
2433 & write (iout,*) 'ran_num',r1
2434 if (r1.lt.0.0d0) OKRandom=.false.
2436 if (.not.OKRandom) then
2437 write (iout,*) 'PRNG IS NOT WORKING!!!'
2438 print *,'PRNG IS NOT WORKING!!!'
2441 call mpi_abort(mpi_comm_world,error_msg,ierr)
2444 write (iout,*) 'too many processors for parallel prng'
2445 write (*,*) 'too many processors for parallel prng'
2453 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)