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 call readi(controlcard,'SYM',symetr,1)
100 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
101 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
102 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
103 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
104 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
105 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
106 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
107 call reada(controlcard,'DRMS',drms,0.1D0)
108 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
109 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
110 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
111 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
112 write (iout,'(a,f10.1)')'DRMS = ',drms
113 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
114 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
116 call readi(controlcard,'NZ_START',nz_start,0)
117 call readi(controlcard,'NZ_END',nz_end,0)
118 call readi(controlcard,'IZ_SC',iz_sc,0)
120 safety = 60.0d0*safety
123 call reada(controlcard,"T_BATH",t_bath,300.0d0)
124 minim=(index(controlcard,'MINIMIZE').gt.0)
125 dccart=(index(controlcard,'CART').gt.0)
126 overlapsc=(index(controlcard,'OVERLAP').gt.0)
127 overlapsc=.not.overlapsc
128 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
129 searchsc=.not.searchsc
130 sideadd=(index(controlcard,'SIDEADD').gt.0)
131 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
132 outpdb=(index(controlcard,'PDBOUT').gt.0)
133 outmol2=(index(controlcard,'MOL2OUT').gt.0)
134 pdbref=(index(controlcard,'PDBREF').gt.0)
135 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
136 indpdb=index(controlcard,'PDBSTART')
137 extconf=(index(controlcard,'EXTCONF').gt.0)
138 call readi(controlcard,'IPRINT',iprint,0)
139 call readi(controlcard,'MAXGEN',maxgen,10000)
140 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
141 call readi(controlcard,"KDIAG",kdiag,0)
142 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
143 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
144 & write (iout,*) "RESCALE_MODE",rescale_mode
145 split_ene=index(controlcard,'SPLIT_ENE').gt.0
146 if (index(controlcard,'REGULAR').gt.0.0D0) then
147 call reada(controlcard,'WEIDIS',weidis,0.1D0)
151 if (index(controlcard,'CHECKGRAD').gt.0) then
153 if (index(controlcard,'CART').gt.0) then
155 elseif (index(controlcard,'CARINT').gt.0) then
160 elseif (index(controlcard,'THREAD').gt.0) then
162 call readi(controlcard,'THREAD',nthread,0)
163 if (nthread.gt.0) then
164 call reada(controlcard,'WEIDIS',weidis,0.1D0)
167 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
168 stop 'Error termination in Read_Control.'
170 else if (index(controlcard,'MCMA').gt.0) then
172 else if (index(controlcard,'MCEE').gt.0) then
174 else if (index(controlcard,'MULTCONF').gt.0) then
176 else if (index(controlcard,'MAP').gt.0) then
178 call readi(controlcard,'MAP',nmap,0)
179 else if (index(controlcard,'CSA').gt.0) then
181 crc else if (index(controlcard,'ZSCORE').gt.0) then
183 crc ZSCORE is rm from UNRES, modecalc=9 is available
186 cfcm else if (index(controlcard,'MCMF').gt.0) then
188 else if (index(controlcard,'SOFTREG').gt.0) then
190 else if (index(controlcard,'CHECK_BOND').gt.0) then
192 else if (index(controlcard,'TEST').gt.0) then
194 else if (index(controlcard,'MD').gt.0) then
196 else if (index(controlcard,'RE ').gt.0) then
200 lmuca=index(controlcard,'MUCA').gt.0
201 call readi(controlcard,'MUCADYN',mucadyn,0)
202 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
203 if (lmuca .and. (me.eq.king .or. .not.out1file ))
205 write (iout,*) 'MUCADYN=',mucadyn
206 write (iout,*) 'MUCASMOOTH=',muca_smooth
209 iscode=index(controlcard,'ONE_LETTER')
210 indphi=index(controlcard,'PHI')
211 indback=index(controlcard,'BACK')
212 iranconf=index(controlcard,'RAND_CONF')
213 i2ndstr=index(controlcard,'USE_SEC_PRED')
214 gradout=index(controlcard,'GRADOUT').gt.0
215 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
216 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
217 if (me.eq.king .or. .not.out1file )
218 & write (iout,*) "DISTCHAINMAX",distchainmax
220 if(me.eq.king.or..not.out1file)
221 & write (iout,'(2a)') diagmeth(kdiag),
222 & ' routine used to diagonalize matrices.'
225 c--------------------------------------------------------------------------
226 subroutine read_REMDpar
230 implicit real*8 (a-h,o-z)
232 include 'COMMON.IOUNITS'
233 include 'COMMON.TIME1'
236 include 'COMMON.LANGEVIN'
238 include 'COMMON.LANGEVIN.lang0'
240 include 'COMMON.INTERACT'
241 include 'COMMON.NAMES'
243 include 'COMMON.REMD'
244 include 'COMMON.CONTROL'
245 include 'COMMON.SETUP'
247 character*320 controlcard
248 character*3200 controlcard1
249 integer iremd_m_total
251 if(me.eq.king.or..not.out1file)
252 & write (iout,*) "REMD setup"
254 call card_concat(controlcard)
255 call readi(controlcard,"NREP",nrep,3)
256 call readi(controlcard,"NSTEX",nstex,1000)
257 call reada(controlcard,"RETMIN",retmin,10.0d0)
258 call reada(controlcard,"RETMAX",retmax,1000.0d0)
259 mremdsync=(index(controlcard,'SYNC').gt.0)
260 call readi(controlcard,"NSYN",i_sync_step,100)
261 restart1file=(index(controlcard,'REST1FILE').gt.0)
262 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
263 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
264 if(max_cache_traj_use.gt.max_cache_traj)
265 & max_cache_traj_use=max_cache_traj
266 if(me.eq.king.or..not.out1file) then
267 cd if (traj1file) then
268 crc caching is in testing - NTWX is not ignored
269 cd write (iout,*) "NTWX value is ignored"
270 cd write (iout,*) " trajectory is stored to one file by master"
271 cd write (iout,*) " before exchange at NSTEX intervals"
273 write (iout,*) "NREP= ",nrep
274 write (iout,*) "NSTEX= ",nstex
275 write (iout,*) "SYNC= ",mremdsync
276 write (iout,*) "NSYN= ",i_sync_step
277 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
280 if (index(controlcard,'TLIST').gt.0) then
282 call card_concat(controlcard1)
283 read(controlcard1,*) (remd_t(i),i=1,nrep)
284 if(me.eq.king.or..not.out1file)
285 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
288 if (index(controlcard,'MLIST').gt.0) then
290 call card_concat(controlcard1)
291 read(controlcard1,*) (remd_m(i),i=1,nrep)
292 if(me.eq.king.or..not.out1file) then
293 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
296 iremd_m_total=iremd_m_total+remd_m(i)
298 write (iout,*) 'Total number of replicas ',iremd_m_total
301 if(me.eq.king.or..not.out1file)
302 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
305 c--------------------------------------------------------------------------
306 subroutine read_MDpar
310 implicit real*8 (a-h,o-z)
312 include 'COMMON.IOUNITS'
313 include 'COMMON.TIME1'
316 include 'COMMON.LANGEVIN'
318 include 'COMMON.LANGEVIN.lang0'
320 include 'COMMON.INTERACT'
321 include 'COMMON.NAMES'
323 include 'COMMON.SETUP'
324 include 'COMMON.CONTROL'
325 include 'COMMON.SPLITELE'
327 character*320 controlcard
329 call card_concat(controlcard)
330 call readi(controlcard,"NSTEP",n_timestep,1000000)
331 call readi(controlcard,"NTWE",ntwe,100)
332 call readi(controlcard,"NTWX",ntwx,1000)
333 call reada(controlcard,"DT",d_time,1.0d-1)
334 call reada(controlcard,"DVMAX",dvmax,2.0d1)
335 call reada(controlcard,"DAMAX",damax,1.0d1)
336 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
337 call readi(controlcard,"LANG",lang,0)
338 RESPA = index(controlcard,"RESPA") .gt. 0
339 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
340 ntime_split0=ntime_split
341 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
342 ntime_split0=ntime_split
343 call reada(controlcard,"R_CUT",r_cut,2.0d0)
344 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
345 rest = index(controlcard,"REST").gt.0
346 tbf = index(controlcard,"TBF").gt.0
347 usampl = index(controlcard,"USAMPL").gt.0
349 mdpdb = index(controlcard,"MDPDB").gt.0
350 call reada(controlcard,"T_BATH",t_bath,300.0d0)
351 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
352 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
353 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
354 if (count_reset_moment.eq.0) count_reset_moment=1000000000
355 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
356 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
357 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
358 if (count_reset_vel.eq.0) count_reset_vel=1000000000
359 large = index(controlcard,"LARGE").gt.0
360 print_compon = index(controlcard,"PRINT_COMPON").gt.0
361 rattle = index(controlcard,"RATTLE").gt.0
362 c if performing umbrella sampling, fragments constrained are read from the fragment file
368 if(me.eq.king.or..not.out1file) then
370 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
372 write (iout,'(a)') "The units are:"
373 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
374 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
375 & " acceleration: angstrom/(48.9 fs)**2"
376 write (iout,'(a)') "energy: kcal/mol, temperature: K"
378 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
379 write (iout,'(a60,f10.5,a)')
380 & "Initial time step of numerical integration:",d_time,
382 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
384 write (iout,'(2a,i4,a)')
385 & "A-MTS algorithm used; initial time step for fast-varying",
386 & " short-range forces split into",ntime_split," steps."
387 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
388 & r_cut," lambda",rlamb
390 write (iout,'(2a,f10.5)')
391 & "Maximum acceleration threshold to reduce the time step",
392 & "/increase split number:",damax
393 write (iout,'(2a,f10.5)')
394 & "Maximum predicted energy drift to reduce the timestep",
395 & "/increase split number:",edriftmax
396 write (iout,'(a60,f10.5)')
397 & "Maximum velocity threshold to reduce velocities:",dvmax
398 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
399 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
400 if (rattle) write (iout,'(a60)')
401 & "Rattle algorithm used to constrain the virtual bonds"
405 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
406 call reada(controlcard,"RWAT",rwat,1.4d0)
407 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
408 surfarea=index(controlcard,"SURFAREA").gt.0
409 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
410 if(me.eq.king.or..not.out1file)then
411 write (iout,'(/a,$)') "Langevin dynamics calculation"
414 & " with direct integration of Langevin equations"
415 else if (lang.eq.2) then
416 write (iout,'(a/)') " with TINKER stochasic MD integrator"
417 else if (lang.eq.3) then
418 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
419 else if (lang.eq.4) then
420 write (iout,'(a/)') " in overdamped mode"
422 write (iout,'(//a,i5)')
423 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
426 write (iout,'(a60,f10.5)') "Temperature:",t_bath
427 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
428 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
429 write (iout,'(a60,f10.5)')
430 & "Scaling factor of the friction forces:",scal_fric
431 if (surfarea) write (iout,'(2a,i10,a)')
432 & "Friction coefficients will be scaled by solvent-accessible",
433 & " surface area every",reset_fricmat," steps."
435 c Calculate friction coefficients and bounds of stochastic forces
436 eta=6*pi*cPoise*etawat
437 if(me.eq.king.or..not.out1file)
438 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
440 gamp=scal_fric*(pstok+rwat)*eta
441 stdfp=dsqrt(2*Rb*t_bath/d_time)
443 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
444 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
446 if(me.eq.king.or..not.out1file)then
447 write (iout,'(/2a/)')
448 & "Radii of site types and friction coefficients and std's of",
449 & " stochastic forces of fully exposed sites"
450 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
452 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
453 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
457 if(me.eq.king.or..not.out1file)then
458 write (iout,'(a)') "Berendsen bath calculation"
459 write (iout,'(a60,f10.5)') "Temperature:",t_bath
460 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
462 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
463 & count_reset_moment," steps"
465 & write (iout,'(a,i10,a)')
466 & "Velocities will be reset at random every",count_reset_vel,
470 if(me.eq.king.or..not.out1file)
471 & write (iout,'(a31)') "Microcanonical mode calculation"
473 if(me.eq.king.or..not.out1file)then
474 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
476 write(iout,*) "MD running with constraints."
477 write(iout,*) "Equilibration time ", eq_time, " mtus."
478 write(iout,*) "Constraining ", nfrag," fragments."
479 write(iout,*) "Length of each fragment, weight and q0:"
481 write (iout,*) "Set of restraints #",iset
483 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
484 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
486 write(iout,*) "constraints between ", npair, "fragments."
487 write(iout,*) "constraint pairs, weights and q0:"
489 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
490 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
492 write(iout,*) "angle constraints within ", nfrag_back,
493 & "backbone fragments."
494 write(iout,*) "fragment, weights:"
496 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
497 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
498 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
501 iset=mod(kolor,nset)+1
504 if(me.eq.king.or..not.out1file)
505 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
508 c------------------------------------------------------------------------------
511 C Read molecular data.
513 implicit real*8 (a-h,o-z)
519 include 'COMMON.IOUNITS'
522 include 'COMMON.INTERACT'
523 include 'COMMON.LOCAL'
524 include 'COMMON.NAMES'
525 include 'COMMON.CHAIN'
526 include 'COMMON.FFIELD'
527 include 'COMMON.SBRIDGE'
528 include 'COMMON.HEADER'
529 include 'COMMON.CONTROL'
530 include 'COMMON.DBASE'
531 include 'COMMON.THREAD'
532 include 'COMMON.CONTACTS'
533 include 'COMMON.TORCNSTR'
534 include 'COMMON.TIME1'
535 include 'COMMON.BOUNDS'
537 include 'COMMON.SETUP'
538 character*4 sequence(maxres)
540 double precision x(maxvar)
541 character*256 pdbfile
542 character*320 weightcard
543 character*80 weightcard_t,ucase
544 dimension itype_pdb(maxres)
545 common /pizda/ itype_pdb
546 logical seq_comp,fail
547 double precision energia(0:n_ene)
553 C Read weights of the subsequent energy terms.
554 call card_concat(weightcard)
555 call reada(weightcard,'WLONG',wlong,1.0D0)
556 call reada(weightcard,'WSC',wsc,wlong)
557 call reada(weightcard,'WSCP',wscp,wlong)
558 call reada(weightcard,'WELEC',welec,1.0D0)
559 call reada(weightcard,'WVDWPP',wvdwpp,welec)
560 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
561 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
562 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
563 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
564 call reada(weightcard,'WTURN3',wturn3,1.0D0)
565 call reada(weightcard,'WTURN4',wturn4,1.0D0)
566 call reada(weightcard,'WTURN6',wturn6,1.0D0)
567 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
568 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
569 call reada(weightcard,'WBOND',wbond,1.0D0)
570 call reada(weightcard,'WTOR',wtor,1.0D0)
571 call reada(weightcard,'WTORD',wtor_d,1.0D0)
572 call reada(weightcard,'WANG',wang,1.0D0)
573 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
574 call reada(weightcard,'SCAL14',scal14,0.4D0)
575 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
576 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
577 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
578 call reada(weightcard,'TEMP0',temp0,300.0d0)
579 if (index(weightcard,'SOFT').gt.0) ipot=6
580 C 12/1/95 Added weight for the multi-body term WCORR
581 call reada(weightcard,'WCORRH',wcorr,1.0D0)
582 if (wcorr4.gt.0.0d0) wcorr=wcorr4
602 if(me.eq.king.or..not.out1file)
603 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
604 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
606 10 format (/'Energy-term weights (unscaled):'//
607 & 'WSCC= ',f10.6,' (SC-SC)'/
608 & 'WSCP= ',f10.6,' (SC-p)'/
609 & 'WELEC= ',f10.6,' (p-p electr)'/
610 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
611 & 'WBOND= ',f10.6,' (stretching)'/
612 & 'WANG= ',f10.6,' (bending)'/
613 & 'WSCLOC= ',f10.6,' (SC local)'/
614 & 'WTOR= ',f10.6,' (torsional)'/
615 & 'WTORD= ',f10.6,' (double torsional)'/
616 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
617 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
618 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
619 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
620 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
621 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
622 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
623 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
624 & 'WTURN6= ',f10.6,' (turns, 6th order)')
625 if(me.eq.king.or..not.out1file)then
626 if (wcorr4.gt.0.0d0) then
627 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
628 & 'between contact pairs of peptide groups'
629 write (iout,'(2(a,f5.3/))')
630 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
631 & 'Range of quenching the correlation terms:',2*delt_corr
632 else if (wcorr.gt.0.0d0) then
633 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
634 & 'between contact pairs of peptide groups'
636 write (iout,'(a,f8.3)')
637 & 'Scaling factor of 1,4 SC-p interactions:',scal14
638 write (iout,'(a,f8.3)')
639 & 'General scaling factor of SC-p interactions:',scalscp
641 r0_corr=cutoff_corr-delt_corr
643 aad(i,1)=scalscp*aad(i,1)
644 aad(i,2)=scalscp*aad(i,2)
645 bad(i,1)=scalscp*bad(i,1)
646 bad(i,2)=scalscp*bad(i,2)
648 call rescale_weights(t_bath)
649 if(me.eq.king.or..not.out1file)
650 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
651 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
653 22 format (/'Energy-term weights (scaled):'//
654 & 'WSCC= ',f10.6,' (SC-SC)'/
655 & 'WSCP= ',f10.6,' (SC-p)'/
656 & 'WELEC= ',f10.6,' (p-p electr)'/
657 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
658 & 'WBOND= ',f10.6,' (stretching)'/
659 & 'WANG= ',f10.6,' (bending)'/
660 & 'WSCLOC= ',f10.6,' (SC local)'/
661 & 'WTOR= ',f10.6,' (torsional)'/
662 & 'WTORD= ',f10.6,' (double torsional)'/
663 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
664 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
665 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
666 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
667 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
668 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
669 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
670 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
671 & 'WTURN6= ',f10.6,' (turns, 6th order)')
672 if(me.eq.king.or..not.out1file)
673 & write (iout,*) "Reference temperature for weights calculation:",
675 call reada(weightcard,"D0CM",d0cm,3.78d0)
676 call reada(weightcard,"AKCM",akcm,15.1d0)
677 call reada(weightcard,"AKTH",akth,11.0d0)
678 call reada(weightcard,"AKCT",akct,12.0d0)
679 call reada(weightcard,"V1SS",v1ss,-1.08d0)
680 call reada(weightcard,"V2SS",v2ss,7.61d0)
681 call reada(weightcard,"V3SS",v3ss,13.7d0)
682 call reada(weightcard,"EBR",ebr,-5.50D0)
683 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
685 dyn_ss_mask(i)=.false.
689 dyn_ssbond_ij(i,j)=1.0d300
692 call reada(weightcard,"HT",Ht,0.0D0)
694 ss_depth=ebr/wsc-0.25*eps(1,1)
695 Ht=Ht/wsc-0.25*eps(1,1)
696 akcm=akcm*wstrain/wsc
697 akth=akth*wstrain/wsc
698 akct=akct*wstrain/wsc
699 v1ss=v1ss*wstrain/wsc
700 v2ss=v2ss*wstrain/wsc
701 v3ss=v3ss*wstrain/wsc
703 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
706 if(me.eq.king.or..not.out1file) then
707 write (iout,*) "Parameters of the SS-bond potential:"
708 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
710 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
711 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
712 write (iout,*)" HT",Ht
713 print *,'indpdb=',indpdb,' pdbref=',pdbref
715 if (indpdb.gt.0 .or. pdbref) then
716 read(inp,'(a)') pdbfile
717 if(me.eq.king.or..not.out1file)
718 & write (iout,'(2a)') 'PDB data will be read from file ',
719 & pdbfile(:ilen(pdbfile))
720 open(ipdbin,file=pdbfile,status='old',err=33)
722 33 write (iout,'(a)') 'Error opening PDB file.'
725 c write (iout,*) 'Begin reading pdb data'
728 c write (iout,*) 'Finished reading pdb data'
730 if(me.eq.king.or..not.out1file)
731 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
732 & ' nstart_sup=',nstart_sup
734 itype_pdb(i)=itype(i)
738 nct=nstart_sup+nsup-1
739 call contact(.false.,ncont_ref,icont_ref,co)
742 if(me.eq.king.or..not.out1file)
743 & write(iout,*)'Adding sidechains'
747 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
750 do while (fail.and.nsi.le.maxsi)
751 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
754 if(fail) write(iout,*)'Adding sidechain failed for res ',
755 & i,' after ',nsi,' trials'
760 if (indpdb.eq.0) then
761 C Read sequence if not taken from the pdb file.
763 c print *,'nres=',nres
764 if (iscode.gt.0) then
765 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
767 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
769 C Convert sequence to numeric code
771 itype(i)=rescode(i,sequence(i),iscode)
773 C Assign initial virtual bond lengths
779 vbld(i+nres)=dsc(iabs(itype(i)))
780 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
781 c write (iout,*) "i",i," itype",itype(i),
782 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
786 c print '(20i4)',(itype(i),i=1,nres)
789 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
791 if (itype(i).eq.ntyp1) then
795 else if (iabs(itype(i+1)).ne.20) then
797 else if (iabs(itype(i)).ne.20) then
804 if(me.eq.king.or..not.out1file)then
805 write (iout,*) "ITEL"
807 write (iout,*) i,itype(i),itel(i)
809 print *,'Call Read_Bridge.'
812 C 8/13/98 Set limits to generating the dihedral angles
817 read (inp,*) ndih_constr
818 if (ndih_constr.gt.0) then
820 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
821 if(me.eq.king.or..not.out1file)then
823 & 'There are',ndih_constr,' constraints on phi angles.'
825 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
829 phi0(i)=deg2rad*phi0(i)
830 drange(i)=deg2rad*drange(i)
832 if(me.eq.king.or..not.out1file)
833 & write (iout,*) 'FTORS',ftors
836 phibound(1,ii) = phi0(i)-drange(i)
837 phibound(2,ii) = phi0(i)+drange(i)
844 write (iout,'(a)') 'Boundaries in phi angle sampling:'
846 write (iout,'(a3,i5,2f10.1)')
847 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
853 cd print *,'NNT=',NNT,' NCT=',NCT
854 if (itype(1).eq.ntyp1) nnt=2
855 if (itype(nres).eq.ntyp1) nct=nct-1
857 if(me.eq.king.or..not.out1file)
858 & write (iout,'(a,i3)') 'nsup=',nsup
860 if (nsup.le.(nct-nnt+1)) then
861 do i=0,nct-nnt+1-nsup
862 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
868 & 'Error - sequences to be superposed do not match.'
871 do i=0,nsup-(nct-nnt+1)
872 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
874 nstart_sup=nstart_sup+i
880 & 'Error - sequences to be superposed do not match.'
883 if (nsup.eq.0) nsup=nct-nnt
884 if (nstart_sup.eq.0) nstart_sup=nnt
885 if (nstart_seq.eq.0) nstart_seq=nnt
886 if(me.eq.king.or..not.out1file)
887 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
888 & ' nstart_seq=',nstart_seq
890 c--- Zscore rms -------
891 if (nz_start.eq.0) nz_start=nnt
892 if (nz_end.eq.0 .and. nsup.gt.0) then
894 else if (nz_end.eq.0) then
897 if(me.eq.king.or..not.out1file)then
898 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
899 write (iout,*) 'IZ_SC=',iz_sc
901 c----------------------
904 if (.not.pdbref) then
905 call read_angles(inp,*38)
907 38 write (iout,'(a)') 'Error reading reference structure.'
909 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
910 stop 'Error reading reference structure'
914 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
924 call contact(.true.,ncont_ref,icont_ref,co)
926 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
928 if (constr_dist.gt.0) call read_dist_constr
929 write (iout,*) "After read_dist_constr nhpb",nhpb
931 if(me.eq.king.or..not.out1file)
932 & write (iout,*) 'Contact order:',co
934 if(me.eq.king.or..not.out1file)
935 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
938 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
940 if(me.eq.king.or..not.out1file)
941 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
942 & icont_ref(1,i),' ',
943 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
947 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
948 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
949 & modecalc.ne.10) then
950 C If input structure hasn't been supplied from the PDB file read or generate
952 if (iranconf.eq.0 .and. .not. extconf) then
953 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
954 & write (iout,'(a)') 'Initial geometry will be read in.'
956 read(inp,'(8f10.5)',end=36,err=36)
957 & ((c(l,k),l=1,3),k=1,nres),
958 & ((c(l,k+nres),l=1,3),k=nnt,nct)
959 write (iout,*) "Exit READ_CART"
960 write (iout,'(8f10.5)')
961 & ((c(l,k),l=1,3),k=1,nres),
962 & ((c(l,k+nres),l=1,3),k=nnt,nct)
963 call int_from_cart1(.true.)
964 write (iout,*) "Finish INT_TO_CART"
967 dc(j,i)=c(j,i+1)-c(j,i)
968 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
972 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
974 dc(j,i+nres)=c(j,i+nres)-c(j,i)
975 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
981 call read_angles(inp,*36)
984 36 write (iout,'(a)') 'Error reading angle file.'
986 call mpi_finalize( MPI_COMM_WORLD,IERR )
988 stop 'Error reading angle file.'
990 else if (extconf) then
991 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
992 & write (iout,'(a)') 'Extended chain initial geometry.'
994 theta(i)=90d0*deg2rad
1000 alph(i)=110d0*deg2rad
1003 omeg(i)=-120d0*deg2rad
1004 if (itype(i).le.0) omeg(i)=-omeg(i)
1007 if(me.eq.king.or..not.out1file)
1008 & write (iout,'(a)') 'Random-generated initial geometry.'
1012 if (me.eq.king .or. fg_rank.eq.0 .and. (
1013 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1017 call gen_rand_conf(itmp,*30)
1019 30 write (iout,*) 'Failed to generate random conformation',
1020 & ', itrial=',itrial
1021 write (*,*) 'Processor:',me,
1022 & ' Failed to generate random conformation',
1032 write (iout,'(a,i3,a)') 'Processor:',me,
1033 & ' error in generating random conformation.'
1034 write (*,'(a,i3,a)') 'Processor:',me,
1035 & ' error in generating random conformation.'
1038 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1044 call gen_rand_conf(itmp,*30)
1046 30 write (iout,*) 'Failed to generate random conformation',
1047 & ', itrial=',itrial
1048 write (*,*) 'Failed to generate random conformation',
1049 & ', itrial=',itrial
1051 write (iout,'(a,i3,a)') 'Processor:',me,
1052 & ' error in generating random conformation.'
1053 write (*,'(a,i3,a)') 'Processor:',me,
1054 & ' error in generating random conformation.'
1059 elseif (modecalc.eq.4) then
1060 read (inp,'(a)') intinname
1061 open (intin,file=intinname,status='old',err=333)
1062 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1063 & write (iout,'(a)') 'intinname',intinname
1064 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1066 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1068 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1070 stop 'Error opening angle file.'
1074 C Generate distance constraints, if the PDB structure is to be regularized.
1075 if (nthread.gt.0) then
1076 call read_threadbase
1079 if (me.eq.king .or. .not. out1file)
1081 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1082 write (iout,'(/a,i3,a)')
1083 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1084 write (iout,'(20i4)') (iss(i),i=1,ns)
1086 write(iout,*)"Running with dynamic disulfide-bond formation"
1088 write (iout,'(/a/)') 'Pre-formed links are:'
1094 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1095 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1101 if (ns.gt.0.and.dyn_ss) then
1105 forcon(i-nss)=forcon(i)
1112 dyn_ss_mask(iss(i))=.true.
1115 if (i2ndstr.gt.0) call secstrp2dihc
1116 c call geom_to_var(nvar,x)
1117 c call etotal(energia(0))
1118 c call enerprint(energia(0))
1119 c call briefout(0,etot)
1121 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1122 cd write (iout,'(a)') 'Variable list:'
1123 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1125 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1126 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1127 & 'Processor',myrank,': end reading molecular data.'
1131 c--------------------------------------------------------------------------
1132 logical function seq_comp(itypea,itypeb,length)
1134 integer length,itypea(length),itypeb(length)
1137 if (itypea(i).ne.itypeb(i)) then
1145 c-----------------------------------------------------------------------------
1146 subroutine read_bridge
1147 C Read information about disulfide bridges.
1148 implicit real*8 (a-h,o-z)
1149 include 'DIMENSIONS'
1153 include 'COMMON.IOUNITS'
1154 include 'COMMON.GEO'
1155 include 'COMMON.VAR'
1156 include 'COMMON.INTERACT'
1157 include 'COMMON.LOCAL'
1158 include 'COMMON.NAMES'
1159 include 'COMMON.CHAIN'
1160 include 'COMMON.FFIELD'
1161 include 'COMMON.SBRIDGE'
1162 include 'COMMON.HEADER'
1163 include 'COMMON.CONTROL'
1164 include 'COMMON.DBASE'
1165 include 'COMMON.THREAD'
1166 include 'COMMON.TIME1'
1167 include 'COMMON.SETUP'
1168 C Read bridging residues.
1169 read (inp,*) ns,(iss(i),i=1,ns)
1171 if(me.eq.king.or..not.out1file)
1172 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1173 C Check whether the specified bridging residues are cystines.
1175 if (itype(iss(i)).ne.1) then
1176 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1177 & 'Do you REALLY think that the residue ',
1178 & restyp(itype(iss(i))),i,
1179 & ' can form a disulfide bridge?!!!'
1180 write (*,'(2a,i3,a)')
1181 & 'Do you REALLY think that the residue ',
1182 & restyp(itype(iss(i))),i,
1183 & ' can form a disulfide bridge?!!!'
1185 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1190 C Read preformed bridges.
1192 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1194 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1197 C Check if the residues involved in bridges are in the specified list of
1198 C bridging residues.
1201 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1202 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1203 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1204 & ' contains residues present in other pairs.'
1205 write (*,'(a,i3,a)') 'Disulfide pair',i,
1206 & ' contains residues present in other pairs.'
1208 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1214 if (ihpb(i).eq.iss(j)) goto 10
1216 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1219 if (jhpb(i).eq.iss(j)) goto 20
1221 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1227 ihpb(i)=ihpb(i)+nres
1228 jhpb(i)=jhpb(i)+nres
1234 c----------------------------------------------------------------------------
1235 subroutine read_x(kanal,*)
1236 implicit real*8 (a-h,o-z)
1237 include 'DIMENSIONS'
1238 include 'COMMON.GEO'
1239 include 'COMMON.VAR'
1240 include 'COMMON.CHAIN'
1241 include 'COMMON.IOUNITS'
1242 include 'COMMON.CONTROL'
1243 include 'COMMON.LOCAL'
1244 include 'COMMON.INTERACT'
1245 c Read coordinates from input
1247 read(kanal,'(8f10.5)',end=10,err=10)
1248 & ((c(l,k),l=1,3),k=1,nres),
1249 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1252 c(j,2*nres)=c(j,nres)
1254 call int_from_cart1(.false.)
1257 dc(j,i)=c(j,i+1)-c(j,i)
1258 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1262 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1264 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1265 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1273 c----------------------------------------------------------------------------
1274 subroutine read_threadbase
1275 implicit real*8 (a-h,o-z)
1276 include 'DIMENSIONS'
1277 include 'COMMON.IOUNITS'
1278 include 'COMMON.GEO'
1279 include 'COMMON.VAR'
1280 include 'COMMON.INTERACT'
1281 include 'COMMON.LOCAL'
1282 include 'COMMON.NAMES'
1283 include 'COMMON.CHAIN'
1284 include 'COMMON.FFIELD'
1285 include 'COMMON.SBRIDGE'
1286 include 'COMMON.HEADER'
1287 include 'COMMON.CONTROL'
1288 include 'COMMON.DBASE'
1289 include 'COMMON.THREAD'
1290 include 'COMMON.TIME1'
1291 C Read pattern database for threading.
1292 read (icbase,*) nseq
1294 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1295 & nres_base(2,i),nres_base(3,i)
1296 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1298 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1299 c & nres_base(2,i),nres_base(3,i)
1300 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1304 if (weidis.eq.0.0D0) weidis=0.1D0
1313 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1314 write (iout,'(a,i5)') 'nexcl: ',nexcl
1315 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1318 c------------------------------------------------------------------------------
1319 subroutine setup_var
1320 implicit real*8 (a-h,o-z)
1321 include 'DIMENSIONS'
1322 include 'COMMON.IOUNITS'
1323 include 'COMMON.GEO'
1324 include 'COMMON.VAR'
1325 include 'COMMON.INTERACT'
1326 include 'COMMON.LOCAL'
1327 include 'COMMON.NAMES'
1328 include 'COMMON.CHAIN'
1329 include 'COMMON.FFIELD'
1330 include 'COMMON.SBRIDGE'
1331 include 'COMMON.HEADER'
1332 include 'COMMON.CONTROL'
1333 include 'COMMON.DBASE'
1334 include 'COMMON.THREAD'
1335 include 'COMMON.TIME1'
1336 C Set up variable list.
1342 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1344 ialph(i,1)=nvar+nside
1348 if (indphi.gt.0) then
1350 else if (indback.gt.0) then
1355 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1358 c----------------------------------------------------------------------------
1359 subroutine gen_dist_constr
1360 C Generate CA distance constraints.
1361 implicit real*8 (a-h,o-z)
1362 include 'DIMENSIONS'
1363 include 'COMMON.IOUNITS'
1364 include 'COMMON.GEO'
1365 include 'COMMON.VAR'
1366 include 'COMMON.INTERACT'
1367 include 'COMMON.LOCAL'
1368 include 'COMMON.NAMES'
1369 include 'COMMON.CHAIN'
1370 include 'COMMON.FFIELD'
1371 include 'COMMON.SBRIDGE'
1372 include 'COMMON.HEADER'
1373 include 'COMMON.CONTROL'
1374 include 'COMMON.DBASE'
1375 include 'COMMON.THREAD'
1376 include 'COMMON.TIME1'
1377 dimension itype_pdb(maxres)
1378 common /pizda/ itype_pdb
1380 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1381 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1382 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1384 do i=nstart_sup,nstart_sup+nsup-1
1385 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1386 cd & ' seq_pdb', restyp(itype_pdb(i))
1387 do j=i+2,nstart_sup+nsup-1
1389 ihpb(nhpb)=i+nstart_seq-nstart_sup
1390 jhpb(nhpb)=j+nstart_seq-nstart_sup
1392 dhpb(nhpb)=dist(i,j)
1395 cd write (iout,'(a)') 'Distance constraints:'
1400 cd if (ii.gt.nres) then
1405 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1406 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1407 cd & dhpb(i),forcon(i)
1411 c----------------------------------------------------------------------------
1413 implicit real*8 (a-h,o-z)
1414 include 'DIMENSIONS'
1415 include 'COMMON.MAP'
1416 include 'COMMON.IOUNITS'
1417 character*3 angid(4) /'THE','PHI','ALP','OME'/
1418 character*80 mapcard,ucase
1420 read (inp,'(a)') mapcard
1421 mapcard=ucase(mapcard)
1422 if (index(mapcard,'PHI').gt.0) then
1424 else if (index(mapcard,'THE').gt.0) then
1426 else if (index(mapcard,'ALP').gt.0) then
1428 else if (index(mapcard,'OME').gt.0) then
1431 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1432 stop 'Error - illegal variable spec in MAP card.'
1434 call readi (mapcard,'RES1',res1(imap),0)
1435 call readi (mapcard,'RES2',res2(imap),0)
1436 if (res1(imap).eq.0) then
1437 res1(imap)=res2(imap)
1438 else if (res2(imap).eq.0) then
1439 res2(imap)=res1(imap)
1441 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1443 & 'Error - illegal definition of variable group in MAP.'
1444 stop 'Error - illegal definition of variable group in MAP.'
1446 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1447 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1448 call readi(mapcard,'NSTEP',nstep(imap),0)
1449 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1451 & 'Illegal boundary and/or step size specification in MAP.'
1452 stop 'Illegal boundary and/or step size specification in MAP.'
1457 c----------------------------------------------------------------------------
1459 implicit real*8 (a-h,o-z)
1460 include 'DIMENSIONS'
1461 include 'COMMON.IOUNITS'
1462 include 'COMMON.GEO'
1463 include 'COMMON.CSA'
1464 include 'COMMON.BANK'
1465 include 'COMMON.CONTROL'
1467 character*620 mcmcard
1468 call card_concat(mcmcard)
1470 call readi(mcmcard,'NCONF',nconf,50)
1471 call readi(mcmcard,'NADD',nadd,0)
1472 call readi(mcmcard,'JSTART',jstart,1)
1473 call readi(mcmcard,'JEND',jend,1)
1474 call readi(mcmcard,'NSTMAX',nstmax,500000)
1475 call readi(mcmcard,'N0',n0,1)
1476 call readi(mcmcard,'N1',n1,6)
1477 call readi(mcmcard,'N2',n2,4)
1478 call readi(mcmcard,'N3',n3,0)
1479 call readi(mcmcard,'N4',n4,0)
1480 call readi(mcmcard,'N5',n5,0)
1481 call readi(mcmcard,'N6',n6,10)
1482 call readi(mcmcard,'N7',n7,0)
1483 call readi(mcmcard,'N8',n8,0)
1484 call readi(mcmcard,'N9',n9,0)
1485 call readi(mcmcard,'N14',n14,0)
1486 call readi(mcmcard,'N15',n15,0)
1487 call readi(mcmcard,'N16',n16,0)
1488 call readi(mcmcard,'N17',n17,0)
1489 call readi(mcmcard,'N18',n18,0)
1491 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1493 call readi(mcmcard,'NDIFF',ndiff,2)
1494 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1495 call readi(mcmcard,'IS1',is1,1)
1496 call readi(mcmcard,'IS2',is2,8)
1497 call readi(mcmcard,'NRAN0',nran0,4)
1498 call readi(mcmcard,'NRAN1',nran1,2)
1499 call readi(mcmcard,'IRR',irr,1)
1500 call readi(mcmcard,'NSEED',nseed,20)
1501 call readi(mcmcard,'NTOTAL',ntotal,10000)
1502 call reada(mcmcard,'CUT1',cut1,2.0d0)
1503 call reada(mcmcard,'CUT2',cut2,5.0d0)
1504 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1505 call readi(mcmcard,'ICMAX',icmax,3)
1506 call readi(mcmcard,'IRESTART',irestart,0)
1507 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1510 call reada(mcmcard,'DELE',dele,20.0d0)
1511 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1512 call readi(mcmcard,'IREF',iref,0)
1513 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1514 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1515 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1516 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1517 write (iout,*) "NCONF_IN",nconf_in
1520 c----------------------------------------------------------------------------
1521 cfmc subroutine mcmfread
1522 cfmc implicit real*8 (a-h,o-z)
1523 cfmc include 'DIMENSIONS'
1524 cfmc include 'COMMON.MCMF'
1525 cfmc include 'COMMON.IOUNITS'
1526 cfmc include 'COMMON.GEO'
1527 cfmc character*80 ucase
1528 cfmc character*620 mcmcard
1529 cfmc call card_concat(mcmcard)
1531 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1532 cfmc write(iout,*)'MAXRANT=',maxrant
1533 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1534 cfmc write(iout,*)'MAXFAM=',maxfam
1535 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1536 cfmc write(iout,*)'NNET1=',nnet1
1537 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1538 cfmc write(iout,*)'NNET2=',nnet2
1539 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1540 cfmc write(iout,*)'NNET3=',nnet3
1541 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1542 cfmc write(iout,*)'ILASTT=',ilastt
1543 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1544 cfmc write(iout,*)'MAXSTR=',maxstr
1545 cfmc maxstr_f=maxstr/maxfam
1546 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1547 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1548 cfmc write(iout,*)'NMCMF=',nmcmf
1549 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1550 cfmc write(iout,*)'IFOCUS=',ifocus
1551 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1552 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1553 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1554 cfmc write(iout,*)'INTPRT=',intprt
1555 cfmc call readi(mcmcard,'IPRT',iprt,100)
1556 cfmc write(iout,*)'IPRT=',iprt
1557 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1558 cfmc write(iout,*)'IMAXTR=',imaxtr
1559 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1560 cfmc write(iout,*)'MAXEVEN=',maxeven
1561 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1562 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1563 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1564 cfmc write(iout,*)'INIMIN=',inimin
1565 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1566 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1567 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1568 cfmc write(iout,*)'NTHREAD=',nthread
1569 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1570 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1571 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1572 cfmc write(iout,*)'MAXPERT=',maxpert
1573 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1574 cfmc write(iout,*)'IRMSD=',irmsd
1575 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1576 cfmc write(iout,*)'DENEMIN=',denemin
1577 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1578 cfmc write(iout,*)'RCUT1S=',rcut1s
1579 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1580 cfmc write(iout,*)'RCUT1E=',rcut1e
1581 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1582 cfmc write(iout,*)'RCUT2S=',rcut2s
1583 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1584 cfmc write(iout,*)'RCUT2E=',rcut2e
1585 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1586 cfmc write(iout,*)'DPERT1=',d_pert1
1587 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1588 cfmc write(iout,*)'DPERT1A=',d_pert1a
1589 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1590 cfmc write(iout,*)'DPERT2=',d_pert2
1591 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1592 cfmc write(iout,*)'DPERT2A=',d_pert2a
1593 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1594 cfmc write(iout,*)'DPERT2B=',d_pert2b
1595 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1596 cfmc write(iout,*)'DPERT2C=',d_pert2c
1597 cfmc d_pert1=deg2rad*d_pert1
1598 cfmc d_pert1a=deg2rad*d_pert1a
1599 cfmc d_pert2=deg2rad*d_pert2
1600 cfmc d_pert2a=deg2rad*d_pert2a
1601 cfmc d_pert2b=deg2rad*d_pert2b
1602 cfmc d_pert2c=deg2rad*d_pert2c
1603 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1604 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1605 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1606 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1607 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1608 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1609 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1610 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1611 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1612 cfmc write(iout,*)'RCUTINI=',rcutini
1613 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1614 cfmc write(iout,*)'GRAT=',grat
1615 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1616 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1620 c----------------------------------------------------------------------------
1622 implicit real*8 (a-h,o-z)
1623 include 'DIMENSIONS'
1624 include 'COMMON.MCM'
1625 include 'COMMON.MCE'
1626 include 'COMMON.IOUNITS'
1628 character*320 mcmcard
1629 call card_concat(mcmcard)
1630 call readi(mcmcard,'MAXACC',maxacc,100)
1631 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1632 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1633 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1634 call readi(mcmcard,'MAXREPM',maxrepm,200)
1635 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1636 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1637 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1638 call reada(mcmcard,'E_UP',e_up,5.0D0)
1639 call reada(mcmcard,'DELTE',delte,0.1D0)
1640 call readi(mcmcard,'NSWEEP',nsweep,5)
1641 call readi(mcmcard,'NSTEPH',nsteph,0)
1642 call readi(mcmcard,'NSTEPC',nstepc,0)
1643 call reada(mcmcard,'TMIN',tmin,298.0D0)
1644 call reada(mcmcard,'TMAX',tmax,298.0D0)
1645 call readi(mcmcard,'NWINDOW',nwindow,0)
1646 call readi(mcmcard,'PRINT_MC',print_mc,0)
1647 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1648 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1649 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1650 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1651 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1652 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1653 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1654 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1655 if (nwindow.gt.0) then
1656 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1658 winlen(i)=winend(i)-winstart(i)+1
1661 if (tmax.lt.tmin) tmax=tmin
1662 if (tmax.eq.tmin) then
1666 if (nstepc.gt.0 .and. nsteph.gt.0) then
1667 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1668 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1670 C Probabilities of different move types
1671 sumpro_type(0)=0.0D0
1672 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1673 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1674 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1675 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1676 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1677 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1678 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1680 print *,'i',i,' sumprotype',sumpro_type(i)
1681 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1682 print *,'i',i,' sumprotype',sumpro_type(i)
1686 c----------------------------------------------------------------------------
1687 subroutine read_minim
1688 implicit real*8 (a-h,o-z)
1689 include 'DIMENSIONS'
1690 include 'COMMON.MINIM'
1691 include 'COMMON.IOUNITS'
1693 character*320 minimcard
1694 call card_concat(minimcard)
1695 call readi(minimcard,'MAXMIN',maxmin,2000)
1696 call readi(minimcard,'MAXFUN',maxfun,5000)
1697 call readi(minimcard,'MINMIN',minmin,maxmin)
1698 call readi(minimcard,'MINFUN',minfun,maxmin)
1699 call reada(minimcard,'TOLF',tolf,1.0D-2)
1700 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1701 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1702 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1703 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1704 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1705 & 'Options in energy minimization:'
1706 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1707 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1708 & 'MinMin:',MinMin,' MinFun:',MinFun,
1709 & ' TolF:',TolF,' RTolF:',RTolF
1712 c----------------------------------------------------------------------------
1713 subroutine read_angles(kanal,*)
1714 implicit real*8 (a-h,o-z)
1715 include 'DIMENSIONS'
1716 include 'COMMON.GEO'
1717 include 'COMMON.VAR'
1718 include 'COMMON.CHAIN'
1719 include 'COMMON.IOUNITS'
1720 include 'COMMON.CONTROL'
1721 c Read angles from input
1723 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1724 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1725 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1726 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1729 c 9/7/01 avoid 180 deg valence angle
1730 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1732 theta(i)=deg2rad*theta(i)
1733 phi(i)=deg2rad*phi(i)
1734 alph(i)=deg2rad*alph(i)
1735 omeg(i)=deg2rad*omeg(i)
1740 c----------------------------------------------------------------------------
1741 subroutine reada(rekord,lancuch,wartosc,default)
1743 character*(*) rekord,lancuch
1744 double precision wartosc,default
1747 iread=index(rekord,lancuch)
1748 if (iread.eq.0) then
1752 iread=iread+ilen(lancuch)+1
1753 read (rekord(iread:),*,err=10,end=10) wartosc
1758 c----------------------------------------------------------------------------
1759 subroutine readi(rekord,lancuch,wartosc,default)
1761 character*(*) rekord,lancuch
1762 integer wartosc,default
1765 iread=index(rekord,lancuch)
1766 if (iread.eq.0) then
1770 iread=iread+ilen(lancuch)+1
1771 read (rekord(iread:),*,err=10,end=10) wartosc
1776 c----------------------------------------------------------------------------
1777 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1780 integer tablica(dim),default
1781 character*(*) rekord,lancuch
1788 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1789 if (iread.eq.0) return
1790 iread=iread+ilen(lancuch)+1
1791 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1794 c----------------------------------------------------------------------------
1795 subroutine multreada(rekord,lancuch,tablica,dim,default)
1798 double precision tablica(dim),default
1799 character*(*) rekord,lancuch
1806 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1807 if (iread.eq.0) return
1808 iread=iread+ilen(lancuch)+1
1809 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1812 c----------------------------------------------------------------------------
1813 subroutine openunits
1814 implicit real*8 (a-h,o-z)
1815 include 'DIMENSIONS'
1818 character*16 form,nodename
1821 include 'COMMON.SETUP'
1822 include 'COMMON.IOUNITS'
1824 include 'COMMON.CONTROL'
1825 integer lenpre,lenpot,ilen,lentmp
1827 character*3 out1file_text,ucase
1830 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1831 call getenv_loc("PREFIX",prefix)
1833 call getenv_loc("POT",pot)
1834 call getenv_loc("DIRTMP",tmpdir)
1835 call getenv_loc("CURDIR",curdir)
1836 call getenv_loc("OUT1FILE",out1file_text)
1837 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1838 out1file_text=ucase(out1file_text)
1839 if (out1file_text(1:1).eq."Y") then
1842 out1file=fg_rank.gt.0
1847 if (lentmp.gt.0) then
1848 write (*,'(80(1h!))')
1849 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1850 write (*,'(80(1h!))')
1851 write (*,*)"All output files will be on node /tmp directory."
1853 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1854 if (me.eq.king) then
1855 write (*,*) "The master node is ",nodename
1856 else if (fg_rank.eq.0) then
1857 write (*,*) "I am the CG slave node ",nodename
1859 write (*,*) "I am the FG slave node ",nodename
1862 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1863 lenpre = lentmp+lenpre+1
1865 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1866 C Get the names and open the input files
1867 #if defined(WINIFL) || defined(WINPGI)
1868 open(1,file=pref_orig(:ilen(pref_orig))//
1869 & '.inp',status='old',readonly,shared)
1870 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1871 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1872 C Get parameter filenames and open the parameter files.
1873 call getenv_loc('BONDPAR',bondname)
1874 open (ibond,file=bondname,status='old',readonly,shared)
1875 call getenv_loc('THETPAR',thetname)
1876 open (ithep,file=thetname,status='old',readonly,shared)
1877 call getenv_loc('ROTPAR',rotname)
1878 open (irotam,file=rotname,status='old',readonly,shared)
1879 call getenv_loc('TORPAR',torname)
1880 open (itorp,file=torname,status='old',readonly,shared)
1881 call getenv_loc('TORDPAR',tordname)
1882 open (itordp,file=tordname,status='old',readonly,shared)
1883 call getenv_loc('FOURIER',fouriername)
1884 open (ifourier,file=fouriername,status='old',readonly,shared)
1885 call getenv_loc('ELEPAR',elename)
1886 open (ielep,file=elename,status='old',readonly,shared)
1887 call getenv_loc('SIDEPAR',sidename)
1888 open (isidep,file=sidename,status='old',readonly,shared)
1889 #elif (defined CRAY) || (defined AIX)
1890 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1892 c print *,"Processor",myrank," opened file 1"
1893 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1894 c print *,"Processor",myrank," opened file 9"
1895 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1896 C Get parameter filenames and open the parameter files.
1897 call getenv_loc('BONDPAR',bondname)
1898 open (ibond,file=bondname,status='old',action='read')
1899 c print *,"Processor",myrank," opened file IBOND"
1900 call getenv_loc('THETPAR',thetname)
1901 open (ithep,file=thetname,status='old',action='read')
1902 c print *,"Processor",myrank," opened file ITHEP"
1903 call getenv_loc('ROTPAR',rotname)
1904 open (irotam,file=rotname,status='old',action='read')
1905 c print *,"Processor",myrank," opened file IROTAM"
1906 call getenv_loc('TORPAR',torname)
1907 open (itorp,file=torname,status='old',action='read')
1908 c print *,"Processor",myrank," opened file ITORP"
1909 call getenv_loc('TORDPAR',tordname)
1910 open (itordp,file=tordname,status='old',action='read')
1911 c print *,"Processor",myrank," opened file ITORDP"
1912 call getenv_loc('SCCORPAR',sccorname)
1913 open (isccor,file=sccorname,status='old',action='read')
1914 c print *,"Processor",myrank," opened file ISCCOR"
1915 call getenv_loc('FOURIER',fouriername)
1916 open (ifourier,file=fouriername,status='old',action='read')
1917 c print *,"Processor",myrank," opened file IFOURIER"
1918 call getenv_loc('ELEPAR',elename)
1919 open (ielep,file=elename,status='old',action='read')
1920 c print *,"Processor",myrank," opened file IELEP"
1921 call getenv_loc('SIDEPAR',sidename)
1922 open (isidep,file=sidename,status='old',action='read')
1923 c print *,"Processor",myrank," opened file ISIDEP"
1924 c print *,"Processor",myrank," opened parameter files"
1926 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1927 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1928 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1929 C Get parameter filenames and open the parameter files.
1930 call getenv_loc('BONDPAR',bondname)
1931 open (ibond,file=bondname,status='old')
1932 call getenv_loc('THETPAR',thetname)
1933 open (ithep,file=thetname,status='old')
1934 call getenv_loc('ROTPAR',rotname)
1935 open (irotam,file=rotname,status='old')
1936 call getenv_loc('TORPAR',torname)
1937 open (itorp,file=torname,status='old')
1938 call getenv_loc('TORDPAR',tordname)
1939 open (itordp,file=tordname,status='old')
1940 call getenv_loc('SCCORPAR',sccorname)
1941 open (isccor,file=sccorname,status='old')
1942 call getenv_loc('FOURIER',fouriername)
1943 open (ifourier,file=fouriername,status='old')
1944 call getenv_loc('ELEPAR',elename)
1945 open (ielep,file=elename,status='old')
1946 call getenv_loc('SIDEPAR',sidename)
1947 open (isidep,file=sidename,status='old')
1949 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1951 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1952 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1953 C Get parameter filenames and open the parameter files.
1954 call getenv_loc('BONDPAR',bondname)
1955 open (ibond,file=bondname,status='old',readonly)
1956 call getenv_loc('THETPAR',thetname)
1957 open (ithep,file=thetname,status='old',readonly)
1958 call getenv_loc('ROTPAR',rotname)
1959 open (irotam,file=rotname,status='old',readonly)
1960 call getenv_loc('TORPAR',torname)
1961 open (itorp,file=torname,status='old',readonly)
1962 call getenv_loc('TORDPAR',tordname)
1963 open (itordp,file=tordname,status='old',readonly)
1964 call getenv_loc('SCCORPAR',sccorname)
1965 open (isccor,file=sccorname,status='old',readonly)
1967 call getenv_loc('THETPARPDB',thetname_pdb)
1968 print *,"thetname_pdb ",thetname_pdb
1969 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1970 print *,ithep_pdb," opened"
1972 call getenv_loc('FOURIER',fouriername)
1973 open (ifourier,file=fouriername,status='old',readonly)
1974 call getenv_loc('ELEPAR',elename)
1975 open (ielep,file=elename,status='old',readonly)
1976 call getenv_loc('SIDEPAR',sidename)
1977 open (isidep,file=sidename,status='old',readonly)
1979 call getenv_loc('ROTPARPDB',rotname_pdb)
1980 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
1985 C 8/9/01 In the newest version SCp interaction constants are read from a file
1986 C Use -DOLDSCP to use hard-coded constants instead.
1988 call getenv_loc('SCPPAR',scpname)
1989 #if defined(WINIFL) || defined(WINPGI)
1990 open (iscpp,file=scpname,status='old',readonly,shared)
1991 #elif (defined CRAY) || (defined AIX)
1992 open (iscpp,file=scpname,status='old',action='read')
1994 open (iscpp,file=scpname,status='old')
1996 open (iscpp,file=scpname,status='old',readonly)
1999 call getenv_loc('PATTERN',patname)
2000 #if defined(WINIFL) || defined(WINPGI)
2001 open (icbase,file=patname,status='old',readonly,shared)
2002 #elif (defined CRAY) || (defined AIX)
2003 open (icbase,file=patname,status='old',action='read')
2005 open (icbase,file=patname,status='old')
2007 open (icbase,file=patname,status='old',readonly)
2010 C Open output file only for CG processes
2011 c print *,"Processor",myrank," fg_rank",fg_rank
2012 if (fg_rank.eq.0) then
2014 if (nodes.eq.1) then
2017 npos = dlog10(dfloat(nodes-1))+1
2019 if (npos.lt.3) npos=3
2020 write (liczba,'(i1)') npos
2021 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2023 write (liczba,form) me
2024 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2025 & liczba(:ilen(liczba))
2026 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2028 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2030 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2031 & liczba(:ilen(liczba))//'.mol2'
2032 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2033 & liczba(:ilen(liczba))//'.stat'
2035 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2036 & //liczba(:ilen(liczba))//'.stat')
2037 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2040 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2041 & liczba(:ilen(liczba))//'.const'
2046 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2047 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2048 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2049 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2050 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2052 & call copy_to_tmp(pref_orig(:ile(pref_orig))//'_'//pot(:lenpot)//
2054 rest2name=prefix(:ilen(prefix))//'.rst'
2056 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2059 #if defined(AIX) || defined(PGI)
2060 if (me.eq.king .or. .not. out1file)
2061 & open(iout,file=outname,status='unknown')
2063 if (fg_rank.gt.0) then
2064 write (liczba,'(i3.3)') myrank/nfgtasks
2065 write (ll,'(bz,i3.3)') fg_rank
2066 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2071 open(igeom,file=intname,status='unknown',position='append')
2072 open(ipdb,file=pdbname,status='unknown')
2073 open(imol2,file=mol2name,status='unknown')
2074 open(istat,file=statname,status='unknown',position='append')
2076 c1out open(iout,file=outname,status='unknown')
2079 if (me.eq.king .or. .not.out1file)
2080 & open(iout,file=outname,status='unknown')
2082 if (fg_rank.gt.0) then
2083 write (liczba,'(i3.3)') myrank/nfgtasks
2084 write (ll,'(bz,i3.3)') fg_rank
2085 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2090 open(igeom,file=intname,status='unknown',access='append')
2091 open(ipdb,file=pdbname,status='unknown')
2092 open(imol2,file=mol2name,status='unknown')
2093 open(istat,file=statname,status='unknown',access='append')
2095 c1out open(iout,file=outname,status='unknown')
2098 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2099 csa_seed=prefix(:lenpre)//'.CSA.seed'
2100 csa_history=prefix(:lenpre)//'.CSA.history'
2101 csa_bank=prefix(:lenpre)//'.CSA.bank'
2102 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2103 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2104 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2105 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2106 csa_int=prefix(:lenpre)//'.int'
2107 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2108 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2109 csa_in=prefix(:lenpre)//'.CSA.in'
2110 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2113 write (iout,'(80(1h-))')
2114 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2115 write (iout,'(80(1h-))')
2116 write (iout,*) "Input file : ",
2117 & pref_orig(:ilen(pref_orig))//'.inp'
2118 write (iout,*) "Output file : ",
2119 & outname(:ilen(outname))
2121 write (iout,*) "Sidechain potential file : ",
2122 & sidename(:ilen(sidename))
2124 write (iout,*) "SCp potential file : ",
2125 & scpname(:ilen(scpname))
2127 write (iout,*) "Electrostatic potential file : ",
2128 & elename(:ilen(elename))
2129 write (iout,*) "Cumulant coefficient file : ",
2130 & fouriername(:ilen(fouriername))
2131 write (iout,*) "Torsional parameter file : ",
2132 & torname(:ilen(torname))
2133 write (iout,*) "Double torsional parameter file : ",
2134 & tordname(:ilen(tordname))
2135 write (iout,*) "SCCOR parameter file : ",
2136 & sccorname(:ilen(sccorname))
2137 write (iout,*) "Bond & inertia constant file : ",
2138 & bondname(:ilen(bondname))
2139 write (iout,*) "Bending parameter file : ",
2140 & thetname(:ilen(thetname))
2141 write (iout,*) "Rotamer parameter file : ",
2142 & rotname(:ilen(rotname))
2143 write (iout,*) "Thetpdb parameter file : ",
2144 & thetname_pdb(:ilen(thetname_pdb))
2145 write (iout,*) "Threading database : ",
2146 & patname(:ilen(patname))
2148 &write (iout,*)" DIRTMP : ",
2150 write (iout,'(80(1h-))')
2154 c----------------------------------------------------------------------------
2155 subroutine card_concat(card)
2156 implicit real*8 (a-h,o-z)
2157 include 'DIMENSIONS'
2158 include 'COMMON.IOUNITS'
2160 character*80 karta,ucase
2162 read (inp,'(a)') karta
2165 do while (karta(80:80).eq.'&')
2166 card=card(:ilen(card)+1)//karta(:79)
2167 read (inp,'(a)') karta
2170 card=card(:ilen(card)+1)//karta
2173 c----------------------------------------------------------------------------------
2175 implicit real*8 (a-h,o-z)
2176 include 'DIMENSIONS'
2177 include 'COMMON.CHAIN'
2178 include 'COMMON.IOUNITS'
2180 open(irest2,file=rest2name,status='unknown')
2181 read(irest2,*) totT,EK,potE,totE,t_bath
2183 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2186 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2189 read (irest2,*) iset
2194 c---------------------------------------------------------------------------------
2195 subroutine read_fragments
2196 implicit real*8 (a-h,o-z)
2197 include 'DIMENSIONS'
2201 include 'COMMON.SETUP'
2202 include 'COMMON.CHAIN'
2203 include 'COMMON.IOUNITS'
2205 include 'COMMON.CONTROL'
2206 read(inp,*) nset,nfrag,npair,nfrag_back
2207 if(me.eq.king.or..not.out1file)
2208 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2209 & " nfrag_back",nfrag_back
2211 read(inp,*) mset(iset)
2213 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2215 if(me.eq.king.or..not.out1file)
2216 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2217 & ifrag(2,i,iset), qinfrag(i,iset)
2220 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2222 if(me.eq.king.or..not.out1file)
2223 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2224 & ipair(2,i,iset), qinpair(i,iset)
2227 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2228 & wfrag_back(3,i,iset),
2229 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2230 if(me.eq.king.or..not.out1file)
2231 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2232 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2237 c-------------------------------------------------------------------------------
2238 subroutine read_dist_constr
2239 implicit real*8 (a-h,o-z)
2240 include 'DIMENSIONS'
2244 include 'COMMON.SETUP'
2245 include 'COMMON.CONTROL'
2246 include 'COMMON.CHAIN'
2247 include 'COMMON.IOUNITS'
2248 include 'COMMON.SBRIDGE'
2249 integer ifrag_(2,100),ipair_(2,100)
2250 double precision wfrag_(100),wpair_(100)
2251 character*500 controlcard
2252 c write (iout,*) "Calling read_dist_constr"
2253 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2255 call card_concat(controlcard)
2256 call readi(controlcard,"NFRAG",nfrag_,0)
2257 call readi(controlcard,"NPAIR",npair_,0)
2258 call readi(controlcard,"NDIST",ndist_,0)
2259 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2260 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2261 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2262 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2263 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2264 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2265 c write (iout,*) "IFRAG"
2267 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2269 c write (iout,*) "IPAIR"
2271 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2275 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2276 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2277 & ifrag_(2,i)=nstart_sup+nsup-1
2278 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2280 if (wfrag_(i).gt.0.0d0) then
2281 do j=ifrag_(1,i),ifrag_(2,i)-1
2282 do k=j+1,ifrag_(2,i)
2283 c write (iout,*) "j",j," k",k
2285 if (constr_dist.eq.1) then
2290 forcon(nhpb)=wfrag_(i)
2291 else if (constr_dist.eq.2) then
2292 if (ddjk.le.dist_cut) then
2297 forcon(nhpb)=wfrag_(i)
2304 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2307 if (.not.out1file .or. me.eq.king)
2308 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2309 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2311 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2312 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2319 if (wpair_(i).gt.0.0d0) then
2327 do j=ifrag_(1,ii),ifrag_(2,ii)
2328 do k=ifrag_(1,jj),ifrag_(2,jj)
2332 forcon(nhpb)=wpair_(i)
2333 dhpb(nhpb)=dist(j,k)
2335 if (.not.out1file .or. me.eq.king)
2336 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2337 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2339 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2340 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2347 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2348 if (forcon(nhpb+1).gt.0.0d0) then
2350 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2352 if (.not.out1file .or. me.eq.king)
2353 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2354 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2356 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2357 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2364 c-------------------------------------------------------------------------------
2366 subroutine flush(iu)
2371 subroutine flush(iu)
2376 c------------------------------------------------------------------------------
2377 subroutine copy_to_tmp(source)
2378 include "DIMENSIONS"
2379 include "COMMON.IOUNITS"
2380 character*(*) source
2381 character* 256 tmpfile
2385 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2386 inquire(file=tmpfile,exist=ex)
2388 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2389 & " to temporary directory..."
2390 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2391 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2395 c------------------------------------------------------------------------------
2396 subroutine move_from_tmp(source)
2397 include "DIMENSIONS"
2398 include "COMMON.IOUNITS"
2399 character*(*) source
2402 write (*,*) "Moving ",source(:ilen(source)),
2403 & " from temporary directory to working directory"
2404 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2405 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2408 c------------------------------------------------------------------------------
2409 subroutine random_init(seed)
2411 C Initialize random number generator
2413 implicit real*8 (a-h,o-z)
2414 include 'DIMENSIONS'
2417 logical OKRandom, prng_restart
2419 integer iseed_array(4)
2421 include 'COMMON.IOUNITS'
2422 include 'COMMON.TIME1'
2423 include 'COMMON.THREAD'
2424 include 'COMMON.SBRIDGE'
2425 include 'COMMON.CONTROL'
2426 include 'COMMON.MCM'
2427 include 'COMMON.MAP'
2428 include 'COMMON.HEADER'
2429 include 'COMMON.CSA'
2430 include 'COMMON.CHAIN'
2431 include 'COMMON.MUCA'
2433 include 'COMMON.FFIELD'
2434 include 'COMMON.SETUP'
2435 iseed=-dint(dabs(seed))
2436 if (iseed.eq.0) then
2437 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2438 & 'Random seed undefined. The program will stop.'
2439 write (*,'(/80(1h*)/20x,a/80(1h*))')
2440 & 'Random seed undefined. The program will stop.'
2442 call mpi_finalize(mpi_comm_world,ierr)
2444 stop 'Bad random seed.'
2447 if (fg_rank.eq.0) then
2451 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2452 OKRandom = prng_restart(me,iseed)
2455 tmp=65536.0d0**(4-i)
2456 iseed_array(i) = dint(seed/tmp)
2457 seed=seed-iseed_array(i)*tmp
2460 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2461 & (iseed_array(i),i=1,4)
2462 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2463 & (iseed_array(i),i=1,4)
2464 OKRandom = prng_restart(me,iseed_array)
2467 c r1 = prng_next(me)
2468 r1=ran_number(0.0D0,1.0D0)
2470 & write (iout,*) 'ran_num',r1
2471 if (r1.lt.0.0d0) OKRandom=.false.
2473 if (.not.OKRandom) then
2474 write (iout,*) 'PRNG IS NOT WORKING!!!'
2475 print *,'PRNG IS NOT WORKING!!!'
2478 call mpi_abort(mpi_comm_world,error_msg,ierr)
2481 write (iout,*) 'too many processors for parallel prng'
2482 write (*,*) 'too many processors for parallel prng'
2490 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)