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 include 'COMMON.SHIELD'
85 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
86 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
88 character*320 controlcard
93 read (INP,'(a)') titel
94 call card_concat(controlcard)
95 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
96 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
97 call reada(controlcard,'SEED',seed,0.0D0)
98 call random_init(seed)
99 C Set up the time limit (caution! The time must be input in minutes!)
100 read_cart=index(controlcard,'READ_CART').gt.0
101 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
102 C this variable with_theta_constr is the variable which allow to read and execute the
103 C constrains on theta angles WITH_THETA_CONSTR is the keyword
104 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
105 write (iout,*) "with_theta_constr ",with_theta_constr
106 call readi(controlcard,'SYM',symetr,1)
107 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
108 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
109 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
110 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
111 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
112 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
113 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
114 call reada(controlcard,'DRMS',drms,0.1D0)
115 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
116 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
117 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
118 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
119 write (iout,'(a,f10.1)')'DRMS = ',drms
120 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
121 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
123 call readi(controlcard,'NZ_START',nz_start,0)
124 call readi(controlcard,'NZ_END',nz_end,0)
125 call readi(controlcard,'IZ_SC',iz_sc,0)
127 safety = 60.0d0*safety
130 call reada(controlcard,"T_BATH",t_bath,300.0d0)
131 minim=(index(controlcard,'MINIMIZE').gt.0)
132 dccart=(index(controlcard,'CART').gt.0)
133 overlapsc=(index(controlcard,'OVERLAP').gt.0)
134 overlapsc=.not.overlapsc
135 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
136 searchsc=.not.searchsc
137 sideadd=(index(controlcard,'SIDEADD').gt.0)
138 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
139 outpdb=(index(controlcard,'PDBOUT').gt.0)
140 outmol2=(index(controlcard,'MOL2OUT').gt.0)
141 pdbref=(index(controlcard,'PDBREF').gt.0)
142 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
143 indpdb=index(controlcard,'PDBSTART')
144 extconf=(index(controlcard,'EXTCONF').gt.0)
145 AFMlog=(index(controlcard,'AFM'))
146 selfguide=(index(controlcard,'SELFGUIDE'))
147 print *,'AFMlog',AFMlog,selfguide,"KUPA"
148 call readi(controlcard,'IPRINT',iprint,0)
149 C SHIELD keyword sets if the shielding effect of side-chains is used
150 C 0 denots no shielding is used all peptide are equally despite the
151 C solvent accesible area
152 C 1 the newly introduced function
153 C 2 reseved for further possible developement
154 call readi(controlcard,'SHIELD',shield_mode,0)
155 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
156 write(iout,*) "shield_mode",shield_mode
158 call readi(controlcard,'TORMODE',tor_mode,0)
159 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
160 write(iout,*) "torsional and valence angle mode",tor_mode
161 call readi(controlcard,'MAXGEN',maxgen,10000)
162 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
163 call readi(controlcard,"KDIAG",kdiag,0)
164 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
165 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
166 & write (iout,*) "RESCALE_MODE",rescale_mode
167 split_ene=index(controlcard,'SPLIT_ENE').gt.0
168 if (index(controlcard,'REGULAR').gt.0.0D0) then
169 call reada(controlcard,'WEIDIS',weidis,0.1D0)
173 if (index(controlcard,'CHECKGRAD').gt.0) then
175 if (index(controlcard,'CART').gt.0) then
177 elseif (index(controlcard,'CARINT').gt.0) then
182 elseif (index(controlcard,'THREAD').gt.0) then
184 call readi(controlcard,'THREAD',nthread,0)
185 if (nthread.gt.0) then
186 call reada(controlcard,'WEIDIS',weidis,0.1D0)
189 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
190 stop 'Error termination in Read_Control.'
192 else if (index(controlcard,'MCMA').gt.0) then
194 else if (index(controlcard,'MCEE').gt.0) then
196 else if (index(controlcard,'MULTCONF').gt.0) then
198 else if (index(controlcard,'MAP').gt.0) then
200 call readi(controlcard,'MAP',nmap,0)
201 else if (index(controlcard,'CSA').gt.0) then
203 crc else if (index(controlcard,'ZSCORE').gt.0) then
205 crc ZSCORE is rm from UNRES, modecalc=9 is available
208 cfcm else if (index(controlcard,'MCMF').gt.0) then
210 else if (index(controlcard,'SOFTREG').gt.0) then
212 else if (index(controlcard,'CHECK_BOND').gt.0) then
214 else if (index(controlcard,'TEST').gt.0) then
216 else if (index(controlcard,'MD').gt.0) then
218 else if (index(controlcard,'RE ').gt.0) then
222 lmuca=index(controlcard,'MUCA').gt.0
223 call readi(controlcard,'MUCADYN',mucadyn,0)
224 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
225 if (lmuca .and. (me.eq.king .or. .not.out1file ))
227 write (iout,*) 'MUCADYN=',mucadyn
228 write (iout,*) 'MUCASMOOTH=',muca_smooth
231 iscode=index(controlcard,'ONE_LETTER')
232 indphi=index(controlcard,'PHI')
233 indback=index(controlcard,'BACK')
234 iranconf=index(controlcard,'RAND_CONF')
235 i2ndstr=index(controlcard,'USE_SEC_PRED')
236 gradout=index(controlcard,'GRADOUT').gt.0
237 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
238 C DISTCHAINMAX become obsolete for periodic boundry condition
239 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
240 C Reading the dimensions of box in x,y,z coordinates
241 call reada(controlcard,'BOXX',boxxsize,100.0d0)
242 call reada(controlcard,'BOXY',boxysize,100.0d0)
243 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
244 c Cutoff range for interactions
245 call reada(controlcard,"R_CUT",r_cut,15.0d0)
246 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
247 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
248 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
249 if (lipthick.gt.0.0d0) then
250 bordliptop=(boxzsize+lipthick)/2.0
251 bordlipbot=bordliptop-lipthick
253 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
254 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
255 buflipbot=bordlipbot+lipbufthick
256 bufliptop=bordliptop-lipbufthick
257 if ((lipbufthick*2.0d0).gt.lipthick)
258 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
260 write(iout,*) "bordliptop=",bordliptop
261 write(iout,*) "bordlipbot=",bordlipbot
262 write(iout,*) "bufliptop=",bufliptop
263 write(iout,*) "buflipbot=",buflipbot
264 write (iout,*) "SHIELD MODE",shield_mode
265 if (shield_mode.gt.0) then
267 C VSolvSphere the volume of solving sphere
269 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
270 C there will be no distinction between proline peptide group and normal peptide
271 C group in case of shielding parameters
272 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
273 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
274 write (iout,*) VSolvSphere,VSolvSphere_div
275 C long axis of side chain
277 long_r_sidechain(i)=vbldsc0(1,i)
278 short_r_sidechain(i)=sigma0(i)
282 if (me.eq.king .or. .not.out1file )
283 & write (iout,*) "DISTCHAINMAX",distchainmax
285 if(me.eq.king.or..not.out1file)
286 & write (iout,'(2a)') diagmeth(kdiag),
287 & ' routine used to diagonalize matrices.'
290 c--------------------------------------------------------------------------
291 subroutine read_REMDpar
295 implicit real*8 (a-h,o-z)
297 include 'COMMON.IOUNITS'
298 include 'COMMON.TIME1'
301 include 'COMMON.LANGEVIN'
303 include 'COMMON.LANGEVIN.lang0'
305 include 'COMMON.INTERACT'
306 include 'COMMON.NAMES'
308 include 'COMMON.REMD'
309 include 'COMMON.CONTROL'
310 include 'COMMON.SETUP'
312 character*320 controlcard
313 character*3200 controlcard1
314 integer iremd_m_total
316 if(me.eq.king.or..not.out1file)
317 & write (iout,*) "REMD setup"
319 call card_concat(controlcard)
320 call readi(controlcard,"NREP",nrep,3)
321 call readi(controlcard,"NSTEX",nstex,1000)
322 call reada(controlcard,"RETMIN",retmin,10.0d0)
323 call reada(controlcard,"RETMAX",retmax,1000.0d0)
324 mremdsync=(index(controlcard,'SYNC').gt.0)
325 call readi(controlcard,"NSYN",i_sync_step,100)
326 restart1file=(index(controlcard,'REST1FILE').gt.0)
327 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
328 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
329 if(max_cache_traj_use.gt.max_cache_traj)
330 & max_cache_traj_use=max_cache_traj
331 if(me.eq.king.or..not.out1file) then
332 cd if (traj1file) then
333 crc caching is in testing - NTWX is not ignored
334 cd write (iout,*) "NTWX value is ignored"
335 cd write (iout,*) " trajectory is stored to one file by master"
336 cd write (iout,*) " before exchange at NSTEX intervals"
338 write (iout,*) "NREP= ",nrep
339 write (iout,*) "NSTEX= ",nstex
340 write (iout,*) "SYNC= ",mremdsync
341 write (iout,*) "NSYN= ",i_sync_step
342 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
345 if (index(controlcard,'TLIST').gt.0) then
347 call card_concat(controlcard1)
348 read(controlcard1,*) (remd_t(i),i=1,nrep)
349 if(me.eq.king.or..not.out1file)
350 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
353 if (index(controlcard,'MLIST').gt.0) then
355 call card_concat(controlcard1)
356 read(controlcard1,*) (remd_m(i),i=1,nrep)
357 if(me.eq.king.or..not.out1file) then
358 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
361 iremd_m_total=iremd_m_total+remd_m(i)
363 write (iout,*) 'Total number of replicas ',iremd_m_total
366 if(me.eq.king.or..not.out1file)
367 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
370 c--------------------------------------------------------------------------
371 subroutine read_MDpar
375 implicit real*8 (a-h,o-z)
377 include 'COMMON.IOUNITS'
378 include 'COMMON.TIME1'
381 include 'COMMON.LANGEVIN'
383 include 'COMMON.LANGEVIN.lang0'
385 include 'COMMON.INTERACT'
386 include 'COMMON.NAMES'
388 include 'COMMON.SETUP'
389 include 'COMMON.CONTROL'
390 include 'COMMON.SPLITELE'
392 character*320 controlcard
394 call card_concat(controlcard)
395 call readi(controlcard,"NSTEP",n_timestep,1000000)
396 call readi(controlcard,"NTWE",ntwe,100)
397 call readi(controlcard,"NTWX",ntwx,1000)
398 call reada(controlcard,"DT",d_time,1.0d-1)
399 call reada(controlcard,"DVMAX",dvmax,2.0d1)
400 call reada(controlcard,"DAMAX",damax,1.0d1)
401 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
402 call readi(controlcard,"LANG",lang,0)
403 RESPA = index(controlcard,"RESPA") .gt. 0
404 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
405 ntime_split0=ntime_split
406 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
407 ntime_split0=ntime_split
408 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
409 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
410 rest = index(controlcard,"REST").gt.0
411 tbf = index(controlcard,"TBF").gt.0
412 usampl = index(controlcard,"USAMPL").gt.0
414 mdpdb = index(controlcard,"MDPDB").gt.0
415 call reada(controlcard,"T_BATH",t_bath,300.0d0)
416 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
417 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
418 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
419 if (count_reset_moment.eq.0) count_reset_moment=1000000000
420 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
421 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
422 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
423 if (count_reset_vel.eq.0) count_reset_vel=1000000000
424 large = index(controlcard,"LARGE").gt.0
425 print_compon = index(controlcard,"PRINT_COMPON").gt.0
426 rattle = index(controlcard,"RATTLE").gt.0
427 c if performing umbrella sampling, fragments constrained are read from the fragment file
433 if(me.eq.king.or..not.out1file) then
435 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
437 write (iout,'(a)') "The units are:"
438 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
439 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
440 & " acceleration: angstrom/(48.9 fs)**2"
441 write (iout,'(a)') "energy: kcal/mol, temperature: K"
443 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
444 write (iout,'(a60,f10.5,a)')
445 & "Initial time step of numerical integration:",d_time,
447 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
449 write (iout,'(2a,i4,a)')
450 & "A-MTS algorithm used; initial time step for fast-varying",
451 & " short-range forces split into",ntime_split," steps."
452 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
453 & r_cut," lambda",rlamb
455 write (iout,'(2a,f10.5)')
456 & "Maximum acceleration threshold to reduce the time step",
457 & "/increase split number:",damax
458 write (iout,'(2a,f10.5)')
459 & "Maximum predicted energy drift to reduce the timestep",
460 & "/increase split number:",edriftmax
461 write (iout,'(a60,f10.5)')
462 & "Maximum velocity threshold to reduce velocities:",dvmax
463 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
464 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
465 if (rattle) write (iout,'(a60)')
466 & "Rattle algorithm used to constrain the virtual bonds"
470 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
471 call reada(controlcard,"RWAT",rwat,1.4d0)
472 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
473 surfarea=index(controlcard,"SURFAREA").gt.0
474 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
475 if(me.eq.king.or..not.out1file)then
476 write (iout,'(/a,$)') "Langevin dynamics calculation"
479 & " with direct integration of Langevin equations"
480 else if (lang.eq.2) then
481 write (iout,'(a/)') " with TINKER stochasic MD integrator"
482 else if (lang.eq.3) then
483 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
484 else if (lang.eq.4) then
485 write (iout,'(a/)') " in overdamped mode"
487 write (iout,'(//a,i5)')
488 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
491 write (iout,'(a60,f10.5)') "Temperature:",t_bath
492 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
493 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
494 write (iout,'(a60,f10.5)')
495 & "Scaling factor of the friction forces:",scal_fric
496 if (surfarea) write (iout,'(2a,i10,a)')
497 & "Friction coefficients will be scaled by solvent-accessible",
498 & " surface area every",reset_fricmat," steps."
500 c Calculate friction coefficients and bounds of stochastic forces
501 eta=6*pi*cPoise*etawat
502 if(me.eq.king.or..not.out1file)
503 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
505 gamp=scal_fric*(pstok+rwat)*eta
506 stdfp=dsqrt(2*Rb*t_bath/d_time)
508 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
509 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
511 if(me.eq.king.or..not.out1file)then
512 write (iout,'(/2a/)')
513 & "Radii of site types and friction coefficients and std's of",
514 & " stochastic forces of fully exposed sites"
515 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
517 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
518 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
522 if(me.eq.king.or..not.out1file)then
523 write (iout,'(a)') "Berendsen bath calculation"
524 write (iout,'(a60,f10.5)') "Temperature:",t_bath
525 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
527 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
528 & count_reset_moment," steps"
530 & write (iout,'(a,i10,a)')
531 & "Velocities will be reset at random every",count_reset_vel,
535 if(me.eq.king.or..not.out1file)
536 & write (iout,'(a31)') "Microcanonical mode calculation"
538 if(me.eq.king.or..not.out1file)then
539 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
541 write(iout,*) "MD running with constraints."
542 write(iout,*) "Equilibration time ", eq_time, " mtus."
543 write(iout,*) "Constraining ", nfrag," fragments."
544 write(iout,*) "Length of each fragment, weight and q0:"
546 write (iout,*) "Set of restraints #",iset
548 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
549 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
551 write(iout,*) "constraints between ", npair, "fragments."
552 write(iout,*) "constraint pairs, weights and q0:"
554 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
555 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
557 write(iout,*) "angle constraints within ", nfrag_back,
558 & "backbone fragments."
559 write(iout,*) "fragment, weights:"
561 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
562 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
563 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
566 iset=mod(kolor,nset)+1
569 if(me.eq.king.or..not.out1file)
570 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
573 c------------------------------------------------------------------------------
576 C Read molecular data.
578 implicit real*8 (a-h,o-z)
584 include 'COMMON.IOUNITS'
587 include 'COMMON.INTERACT'
588 include 'COMMON.LOCAL'
589 include 'COMMON.NAMES'
590 include 'COMMON.CHAIN'
591 include 'COMMON.FFIELD'
592 include 'COMMON.SBRIDGE'
593 include 'COMMON.HEADER'
594 include 'COMMON.CONTROL'
595 include 'COMMON.DBASE'
596 include 'COMMON.THREAD'
597 include 'COMMON.CONTACTS'
598 include 'COMMON.TORCNSTR'
599 include 'COMMON.TIME1'
600 include 'COMMON.BOUNDS'
602 include 'COMMON.SETUP'
603 include 'COMMON.SHIELD'
604 character*4 sequence(maxres)
606 double precision x(maxvar)
607 character*256 pdbfile
608 character*400 weightcard
609 character*80 weightcard_t,ucase
610 dimension itype_pdb(maxres)
611 common /pizda/ itype_pdb
612 logical seq_comp,fail
613 double precision energia(0:n_ene)
619 C Read weights of the subsequent energy terms.
620 call card_concat(weightcard)
621 call reada(weightcard,'WLONG',wlong,1.0D0)
622 call reada(weightcard,'WSC',wsc,wlong)
623 call reada(weightcard,'WSCP',wscp,wlong)
624 call reada(weightcard,'WELEC',welec,1.0D0)
625 call reada(weightcard,'WVDWPP',wvdwpp,welec)
626 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
627 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
628 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
629 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
630 call reada(weightcard,'WTURN3',wturn3,1.0D0)
631 call reada(weightcard,'WTURN4',wturn4,1.0D0)
632 call reada(weightcard,'WTURN6',wturn6,1.0D0)
633 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
634 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
635 call reada(weightcard,'WBOND',wbond,1.0D0)
636 call reada(weightcard,'WTOR',wtor,1.0D0)
637 call reada(weightcard,'WTORD',wtor_d,1.0D0)
638 call reada(weightcard,'WANG',wang,1.0D0)
639 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
640 call reada(weightcard,'SCAL14',scal14,0.4D0)
641 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
642 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
643 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
644 call reada(weightcard,'TEMP0',temp0,300.0d0)
645 call reada(weightcard,'WSHIELD',wshield,1.0d0)
646 call reada(weightcard,'WLT',wliptran,0.0D0)
647 if (index(weightcard,'SOFT').gt.0) ipot=6
648 C 12/1/95 Added weight for the multi-body term WCORR
649 call reada(weightcard,'WCORRH',wcorr,1.0D0)
650 if (wcorr4.gt.0.0d0) wcorr=wcorr4
670 if(me.eq.king.or..not.out1file)
671 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
672 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
674 10 format (/'Energy-term weights (unscaled):'//
675 & 'WSCC= ',f10.6,' (SC-SC)'/
676 & 'WSCP= ',f10.6,' (SC-p)'/
677 & 'WELEC= ',f10.6,' (p-p electr)'/
678 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
679 & 'WBOND= ',f10.6,' (stretching)'/
680 & 'WANG= ',f10.6,' (bending)'/
681 & 'WSCLOC= ',f10.6,' (SC local)'/
682 & 'WTOR= ',f10.6,' (torsional)'/
683 & 'WTORD= ',f10.6,' (double torsional)'/
684 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
685 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
686 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
687 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
688 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
689 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
690 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
691 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
692 & 'WTURN6= ',f10.6,' (turns, 6th order)')
693 if(me.eq.king.or..not.out1file)then
694 if (wcorr4.gt.0.0d0) then
695 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
696 & 'between contact pairs of peptide groups'
697 write (iout,'(2(a,f5.3/))')
698 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
699 & 'Range of quenching the correlation terms:',2*delt_corr
700 else if (wcorr.gt.0.0d0) then
701 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
702 & 'between contact pairs of peptide groups'
704 write (iout,'(a,f8.3)')
705 & 'Scaling factor of 1,4 SC-p interactions:',scal14
706 write (iout,'(a,f8.3)')
707 & 'General scaling factor of SC-p interactions:',scalscp
709 r0_corr=cutoff_corr-delt_corr
711 aad(i,1)=scalscp*aad(i,1)
712 aad(i,2)=scalscp*aad(i,2)
713 bad(i,1)=scalscp*bad(i,1)
714 bad(i,2)=scalscp*bad(i,2)
716 call rescale_weights(t_bath)
717 if(me.eq.king.or..not.out1file)
718 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
719 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
721 22 format (/'Energy-term weights (scaled):'//
722 & 'WSCC= ',f10.6,' (SC-SC)'/
723 & 'WSCP= ',f10.6,' (SC-p)'/
724 & 'WELEC= ',f10.6,' (p-p electr)'/
725 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
726 & 'WBOND= ',f10.6,' (stretching)'/
727 & 'WANG= ',f10.6,' (bending)'/
728 & 'WSCLOC= ',f10.6,' (SC local)'/
729 & 'WTOR= ',f10.6,' (torsional)'/
730 & 'WTORD= ',f10.6,' (double torsional)'/
731 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
732 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
733 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
734 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
735 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
736 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
737 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
738 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
739 & 'WTURN6= ',f10.6,' (turns, 6th order)')
740 if(me.eq.king.or..not.out1file)
741 & write (iout,*) "Reference temperature for weights calculation:",
743 call reada(weightcard,"D0CM",d0cm,3.78d0)
744 call reada(weightcard,"AKCM",akcm,15.1d0)
745 call reada(weightcard,"AKTH",akth,11.0d0)
746 call reada(weightcard,"AKCT",akct,12.0d0)
747 call reada(weightcard,"V1SS",v1ss,-1.08d0)
748 call reada(weightcard,"V2SS",v2ss,7.61d0)
749 call reada(weightcard,"V3SS",v3ss,13.7d0)
750 call reada(weightcard,"EBR",ebr,-5.50D0)
751 call reada(weightcard,"ATRISS",atriss,0.301D0)
752 call reada(weightcard,"BTRISS",btriss,0.021D0)
753 call reada(weightcard,"CTRISS",ctriss,1.001D0)
754 call reada(weightcard,"DTRISS",dtriss,1.001D0)
755 write (iout,*) "ATRISS=", atriss
756 write (iout,*) "BTRISS=", btriss
757 write (iout,*) "CTRISS=", ctriss
758 write (iout,*) "DTRISS=", dtriss
759 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
761 dyn_ss_mask(i)=.false.
765 dyn_ssbond_ij(i,j)=1.0d300
768 call reada(weightcard,"HT",Ht,0.0D0)
770 ss_depth=ebr/wsc-0.25*eps(1,1)
771 Ht=Ht/wsc-0.25*eps(1,1)
772 akcm=akcm*wstrain/wsc
773 akth=akth*wstrain/wsc
774 akct=akct*wstrain/wsc
775 v1ss=v1ss*wstrain/wsc
776 v2ss=v2ss*wstrain/wsc
777 v3ss=v3ss*wstrain/wsc
779 if (wstrain.ne.0.0) then
780 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
786 if(me.eq.king.or..not.out1file) then
787 write (iout,*) "Parameters of the SS-bond potential:"
788 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
790 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
791 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
792 write (iout,*)" HT",Ht
793 print *,'indpdb=',indpdb,' pdbref=',pdbref
795 if (indpdb.gt.0 .or. pdbref) then
796 read(inp,'(a)') pdbfile
797 if(me.eq.king.or..not.out1file)
798 & write (iout,'(2a)') 'PDB data will be read from file ',
799 & pdbfile(:ilen(pdbfile))
800 open(ipdbin,file=pdbfile,status='old',err=33)
802 33 write (iout,'(a)') 'Error opening PDB file.'
805 c write (iout,*) 'Begin reading pdb data'
808 c write (iout,*) 'Finished reading pdb data'
810 if(me.eq.king.or..not.out1file)
811 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
812 & ' nstart_sup=',nstart_sup
814 itype_pdb(i)=itype(i)
818 nct=nstart_sup+nsup-1
819 call contact(.false.,ncont_ref,icont_ref,co)
822 if(me.eq.king.or..not.out1file)
823 & write(iout,*)'Adding sidechains'
827 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
830 do while (fail.and.nsi.le.maxsi)
831 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
834 if(fail) write(iout,*)'Adding sidechain failed for res ',
835 & i,' after ',nsi,' trials'
840 if (indpdb.eq.0) then
841 C Read sequence if not taken from the pdb file.
843 c print *,'nres=',nres
844 if (iscode.gt.0) then
845 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
847 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
849 C Convert sequence to numeric code
851 itype(i)=rescode(i,sequence(i),iscode)
853 C Assign initial virtual bond lengths
859 vbld(i+nres)=dsc(iabs(itype(i)))
860 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
861 c write (iout,*) "i",i," itype",itype(i),
862 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
866 c print '(20i4)',(itype(i),i=1,nres)
869 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
871 if (itype(i).eq.ntyp1) then
875 else if (iabs(itype(i+1)).ne.20) then
877 else if (iabs(itype(i)).ne.20) then
884 if(me.eq.king.or..not.out1file)then
885 write (iout,*) "ITEL"
887 write (iout,*) i,itype(i),itel(i)
889 print *,'Call Read_Bridge.'
892 C 8/13/98 Set limits to generating the dihedral angles
897 read (inp,*) ndih_constr
898 if (ndih_constr.gt.0) then
900 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
902 if(me.eq.king.or..not.out1file)then
904 & 'There are',ndih_constr,' constraints on phi angles.'
906 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
911 phi0(i)=deg2rad*phi0(i)
912 drange(i)=deg2rad*drange(i)
914 C if(me.eq.king.or..not.out1file)
915 C & write (iout,*) 'FTORS',ftors
918 phibound(1,ii) = phi0(i)-drange(i)
919 phibound(2,ii) = phi0(i)+drange(i)
922 C first setting the theta boundaries to 0 to pi
923 C this mean that there is no energy penalty for any angle occuring this can be applied
924 C for generate random conformation but is not implemented in this way
929 C begin reading theta constrains this is quartic constrains allowing to
930 C have smooth second derivative
931 if (with_theta_constr) then
932 C with_theta_constr is keyword allowing for occurance of theta constrains
933 read (inp,*) ntheta_constr
934 C ntheta_constr is the number of theta constrains
935 if (ntheta_constr.gt.0) then
937 read (inp,*) (itheta_constr(i),theta_constr0(i),
938 & theta_drange(i),for_thet_constr(i),
940 C the above code reads from 1 to ntheta_constr
941 C itheta_constr(i) residue i for which is theta_constr
942 C theta_constr0 the global minimum value
943 C theta_drange is range for which there is no energy penalty
944 C for_thet_constr is the force constant for quartic energy penalty
946 if(me.eq.king.or..not.out1file)then
948 & 'There are',ntheta_constr,' constraints on phi angles.'
950 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
956 theta_constr0(i)=deg2rad*theta_constr0(i)
957 theta_drange(i)=deg2rad*theta_drange(i)
959 C if(me.eq.king.or..not.out1file)
960 C & write (iout,*) 'FTORS',ftors
961 C do i=1,ntheta_constr
962 C ii = itheta_constr(i)
963 C thetabound(1,ii) = phi0(i)-drange(i)
964 C thetabound(2,ii) = phi0(i)+drange(i)
966 endif ! ntheta_constr.gt.0
967 endif! with_theta_constr
969 C with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
970 C write (iout,*) "with_dihed_constr ",with_dihed_constr
975 write (iout,'(a)') 'Boundaries in phi angle sampling:'
977 write (iout,'(a3,i5,2f10.1)')
978 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
984 cd print *,'NNT=',NNT,' NCT=',NCT
985 if (itype(1).eq.ntyp1) nnt=2
986 if (itype(nres).eq.ntyp1) nct=nct-1
988 if(me.eq.king.or..not.out1file)
989 & write (iout,'(a,i3)') 'nsup=',nsup
991 if (nsup.le.(nct-nnt+1)) then
992 do i=0,nct-nnt+1-nsup
993 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
999 & 'Error - sequences to be superposed do not match.'
1002 do i=0,nsup-(nct-nnt+1)
1003 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1005 nstart_sup=nstart_sup+i
1011 & 'Error - sequences to be superposed do not match.'
1014 if (nsup.eq.0) nsup=nct-nnt
1015 if (nstart_sup.eq.0) nstart_sup=nnt
1016 if (nstart_seq.eq.0) nstart_seq=nnt
1017 if(me.eq.king.or..not.out1file)
1018 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1019 & ' nstart_seq=',nstart_seq
1021 c--- Zscore rms -------
1022 if (nz_start.eq.0) nz_start=nnt
1023 if (nz_end.eq.0 .and. nsup.gt.0) then
1025 else if (nz_end.eq.0) then
1028 if(me.eq.king.or..not.out1file)then
1029 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1030 write (iout,*) 'IZ_SC=',iz_sc
1032 c----------------------
1035 if (.not.pdbref) then
1036 call read_angles(inp,*38)
1038 38 write (iout,'(a)') 'Error reading reference structure.'
1040 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1041 stop 'Error reading reference structure'
1045 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1052 cref(j,i,kkk)=c(j,i)
1055 call contact(.true.,ncont_ref,icont_ref,co)
1059 write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1061 if (constr_dist.gt.0) call read_dist_constr
1062 write (iout,*) "After read_dist_constr nhpb",nhpb
1063 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1065 if(me.eq.king.or..not.out1file)
1066 & write (iout,*) 'Contact order:',co
1068 if(me.eq.king.or..not.out1file)
1069 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1072 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1074 if(me.eq.king.or..not.out1file)
1075 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1076 & icont_ref(1,i),' ',
1077 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1081 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1082 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1083 & modecalc.ne.10) then
1084 C If input structure hasn't been supplied from the PDB file read or generate
1086 if (iranconf.eq.0 .and. .not. extconf) then
1087 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1088 & write (iout,'(a)') 'Initial geometry will be read in.'
1090 read(inp,'(8f10.5)',end=36,err=36)
1091 & ((c(l,k),l=1,3),k=1,nres),
1092 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1093 write (iout,*) "Exit READ_CART"
1094 write (iout,'(8f10.5)')
1095 & ((c(l,k),l=1,3),k=1,nres),
1096 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1097 call int_from_cart1(.true.)
1098 write (iout,*) "Finish INT_TO_CART"
1101 dc(j,i)=c(j,i+1)-c(j,i)
1102 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1106 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1108 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1109 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1115 call read_angles(inp,*36)
1118 36 write (iout,'(a)') 'Error reading angle file.'
1120 call mpi_finalize( MPI_COMM_WORLD,IERR )
1122 stop 'Error reading angle file.'
1124 else if (extconf) then
1125 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1126 & write (iout,'(a)') 'Extended chain initial geometry.'
1128 theta(i)=90d0*deg2rad
1131 phi(i)=180d0*deg2rad
1134 alph(i)=110d0*deg2rad
1137 omeg(i)=-120d0*deg2rad
1138 if (itype(i).le.0) omeg(i)=-omeg(i)
1140 call chainbuild_extconf
1142 if(me.eq.king.or..not.out1file)
1143 & write (iout,'(a)') 'Random-generated initial geometry.'
1147 if (me.eq.king .or. fg_rank.eq.0 .and. (
1148 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1152 call gen_rand_conf(itmp,*30)
1154 30 write (iout,*) 'Failed to generate random conformation',
1155 & ', itrial=',itrial
1156 write (*,*) 'Processor:',me,
1157 & ' Failed to generate random conformation',
1167 write (iout,'(a,i3,a)') 'Processor:',me,
1168 & ' error in generating random conformation.'
1169 write (*,'(a,i3,a)') 'Processor:',me,
1170 & ' error in generating random conformation.'
1173 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1178 & ' error in generating random conformation.'
1183 elseif (modecalc.eq.4) then
1184 read (inp,'(a)') intinname
1185 open (intin,file=intinname,status='old',err=333)
1186 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1187 & write (iout,'(a)') 'intinname',intinname
1188 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1190 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1192 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1194 stop 'Error opening angle file.'
1198 C Generate distance constraints, if the PDB structure is to be regularized.
1199 if (nthread.gt.0) then
1200 call read_threadbase
1203 if (me.eq.king .or. .not. out1file)
1205 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1206 write (iout,'(/a,i3,a)')
1207 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1208 write (iout,'(20i4)') (iss(i),i=1,ns)
1210 write(iout,*)"Running with dynamic disulfide-bond formation"
1212 write (iout,'(/a/)') 'Pre-formed links are:'
1218 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1219 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1225 if (ns.gt.0.and.dyn_ss) then
1229 forcon(i-nss)=forcon(i)
1236 dyn_ss_mask(iss(i))=.true.
1239 if (i2ndstr.gt.0) call secstrp2dihc
1240 c call geom_to_var(nvar,x)
1241 c call etotal(energia(0))
1242 c call enerprint(energia(0))
1243 c call briefout(0,etot)
1245 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1246 cd write (iout,'(a)') 'Variable list:'
1247 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1249 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1250 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1251 & 'Processor',myrank,': end reading molecular data.'
1256 c--------------------------------------------------------------------------
1257 logical function seq_comp(itypea,itypeb,length)
1259 integer length,itypea(length),itypeb(length)
1262 if (itypea(i).ne.itypeb(i)) then
1270 c-----------------------------------------------------------------------------
1271 subroutine read_bridge
1272 C Read information about disulfide bridges.
1273 implicit real*8 (a-h,o-z)
1274 include 'DIMENSIONS'
1278 include 'COMMON.IOUNITS'
1279 include 'COMMON.GEO'
1280 include 'COMMON.VAR'
1281 include 'COMMON.INTERACT'
1282 include 'COMMON.LOCAL'
1283 include 'COMMON.NAMES'
1284 include 'COMMON.CHAIN'
1285 include 'COMMON.FFIELD'
1286 include 'COMMON.SBRIDGE'
1287 include 'COMMON.HEADER'
1288 include 'COMMON.CONTROL'
1289 include 'COMMON.DBASE'
1290 include 'COMMON.THREAD'
1291 include 'COMMON.TIME1'
1292 include 'COMMON.SETUP'
1293 C Read bridging residues.
1294 read (inp,*) ns,(iss(i),i=1,ns)
1296 if(me.eq.king.or..not.out1file)
1297 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1298 C Check whether the specified bridging residues are cystines.
1300 if (itype(iss(i)).ne.1) then
1301 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1302 & 'Do you REALLY think that the residue ',
1303 & restyp(itype(iss(i))),i,
1304 & ' can form a disulfide bridge?!!!'
1305 write (*,'(2a,i3,a)')
1306 & 'Do you REALLY think that the residue ',
1307 & restyp(itype(iss(i))),i,
1308 & ' can form a disulfide bridge?!!!'
1310 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1315 C Read preformed bridges.
1317 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1319 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1322 C Check if the residues involved in bridges are in the specified list of
1323 C bridging residues.
1326 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1327 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1328 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1329 & ' contains residues present in other pairs.'
1330 write (*,'(a,i3,a)') 'Disulfide pair',i,
1331 & ' contains residues present in other pairs.'
1333 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1339 if (ihpb(i).eq.iss(j)) goto 10
1341 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1344 if (jhpb(i).eq.iss(j)) goto 20
1346 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1352 ihpb(i)=ihpb(i)+nres
1353 jhpb(i)=jhpb(i)+nres
1359 c----------------------------------------------------------------------------
1360 subroutine read_x(kanal,*)
1361 implicit real*8 (a-h,o-z)
1362 include 'DIMENSIONS'
1363 include 'COMMON.GEO'
1364 include 'COMMON.VAR'
1365 include 'COMMON.CHAIN'
1366 include 'COMMON.IOUNITS'
1367 include 'COMMON.CONTROL'
1368 include 'COMMON.LOCAL'
1369 include 'COMMON.INTERACT'
1370 c Read coordinates from input
1372 read(kanal,'(8f10.5)',end=10,err=10)
1373 & ((c(l,k),l=1,3),k=1,nres),
1374 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1377 c(j,2*nres)=c(j,nres)
1379 call int_from_cart1(.false.)
1382 dc(j,i)=c(j,i+1)-c(j,i)
1383 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1387 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1389 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1390 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1398 c----------------------------------------------------------------------------
1399 subroutine read_threadbase
1400 implicit real*8 (a-h,o-z)
1401 include 'DIMENSIONS'
1402 include 'COMMON.IOUNITS'
1403 include 'COMMON.GEO'
1404 include 'COMMON.VAR'
1405 include 'COMMON.INTERACT'
1406 include 'COMMON.LOCAL'
1407 include 'COMMON.NAMES'
1408 include 'COMMON.CHAIN'
1409 include 'COMMON.FFIELD'
1410 include 'COMMON.SBRIDGE'
1411 include 'COMMON.HEADER'
1412 include 'COMMON.CONTROL'
1413 include 'COMMON.DBASE'
1414 include 'COMMON.THREAD'
1415 include 'COMMON.TIME1'
1416 C Read pattern database for threading.
1417 read (icbase,*) nseq
1419 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1420 & nres_base(2,i),nres_base(3,i)
1421 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1423 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1424 c & nres_base(2,i),nres_base(3,i)
1425 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1429 if (weidis.eq.0.0D0) weidis=0.1D0
1438 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1439 write (iout,'(a,i5)') 'nexcl: ',nexcl
1440 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1443 c------------------------------------------------------------------------------
1444 subroutine setup_var
1445 implicit real*8 (a-h,o-z)
1446 include 'DIMENSIONS'
1447 include 'COMMON.IOUNITS'
1448 include 'COMMON.GEO'
1449 include 'COMMON.VAR'
1450 include 'COMMON.INTERACT'
1451 include 'COMMON.LOCAL'
1452 include 'COMMON.NAMES'
1453 include 'COMMON.CHAIN'
1454 include 'COMMON.FFIELD'
1455 include 'COMMON.SBRIDGE'
1456 include 'COMMON.HEADER'
1457 include 'COMMON.CONTROL'
1458 include 'COMMON.DBASE'
1459 include 'COMMON.THREAD'
1460 include 'COMMON.TIME1'
1461 C Set up variable list.
1467 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1469 ialph(i,1)=nvar+nside
1473 if (indphi.gt.0) then
1475 else if (indback.gt.0) then
1480 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1483 c----------------------------------------------------------------------------
1484 subroutine gen_dist_constr
1485 C Generate CA distance constraints.
1486 implicit real*8 (a-h,o-z)
1487 include 'DIMENSIONS'
1488 include 'COMMON.IOUNITS'
1489 include 'COMMON.GEO'
1490 include 'COMMON.VAR'
1491 include 'COMMON.INTERACT'
1492 include 'COMMON.LOCAL'
1493 include 'COMMON.NAMES'
1494 include 'COMMON.CHAIN'
1495 include 'COMMON.FFIELD'
1496 include 'COMMON.SBRIDGE'
1497 include 'COMMON.HEADER'
1498 include 'COMMON.CONTROL'
1499 include 'COMMON.DBASE'
1500 include 'COMMON.THREAD'
1501 include 'COMMON.TIME1'
1502 dimension itype_pdb(maxres)
1503 common /pizda/ itype_pdb
1505 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1506 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1507 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1509 do i=nstart_sup,nstart_sup+nsup-1
1510 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1511 cd & ' seq_pdb', restyp(itype_pdb(i))
1512 do j=i+2,nstart_sup+nsup-1
1514 ihpb(nhpb)=i+nstart_seq-nstart_sup
1515 jhpb(nhpb)=j+nstart_seq-nstart_sup
1517 dhpb(nhpb)=dist(i,j)
1520 cd write (iout,'(a)') 'Distance constraints:'
1525 cd if (ii.gt.nres) then
1530 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1531 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1532 cd & dhpb(i),forcon(i)
1536 c----------------------------------------------------------------------------
1538 implicit real*8 (a-h,o-z)
1539 include 'DIMENSIONS'
1540 include 'COMMON.MAP'
1541 include 'COMMON.IOUNITS'
1542 character*3 angid(4) /'THE','PHI','ALP','OME'/
1543 character*80 mapcard,ucase
1545 read (inp,'(a)') mapcard
1546 mapcard=ucase(mapcard)
1547 if (index(mapcard,'PHI').gt.0) then
1549 else if (index(mapcard,'THE').gt.0) then
1551 else if (index(mapcard,'ALP').gt.0) then
1553 else if (index(mapcard,'OME').gt.0) then
1556 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1557 stop 'Error - illegal variable spec in MAP card.'
1559 call readi (mapcard,'RES1',res1(imap),0)
1560 call readi (mapcard,'RES2',res2(imap),0)
1561 if (res1(imap).eq.0) then
1562 res1(imap)=res2(imap)
1563 else if (res2(imap).eq.0) then
1564 res2(imap)=res1(imap)
1566 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1568 & 'Error - illegal definition of variable group in MAP.'
1569 stop 'Error - illegal definition of variable group in MAP.'
1571 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1572 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1573 call readi(mapcard,'NSTEP',nstep(imap),0)
1574 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1576 & 'Illegal boundary and/or step size specification in MAP.'
1577 stop 'Illegal boundary and/or step size specification in MAP.'
1582 c----------------------------------------------------------------------------
1584 implicit real*8 (a-h,o-z)
1585 include 'DIMENSIONS'
1586 include 'COMMON.IOUNITS'
1587 include 'COMMON.GEO'
1588 include 'COMMON.CSA'
1589 include 'COMMON.BANK'
1590 include 'COMMON.CONTROL'
1592 character*620 mcmcard
1593 call card_concat(mcmcard)
1595 call readi(mcmcard,'NCONF',nconf,50)
1596 call readi(mcmcard,'NADD',nadd,0)
1597 call readi(mcmcard,'JSTART',jstart,1)
1598 call readi(mcmcard,'JEND',jend,1)
1599 call readi(mcmcard,'NSTMAX',nstmax,500000)
1600 call readi(mcmcard,'N0',n0,1)
1601 call readi(mcmcard,'N1',n1,6)
1602 call readi(mcmcard,'N2',n2,4)
1603 call readi(mcmcard,'N3',n3,0)
1604 call readi(mcmcard,'N4',n4,0)
1605 call readi(mcmcard,'N5',n5,0)
1606 call readi(mcmcard,'N6',n6,10)
1607 call readi(mcmcard,'N7',n7,0)
1608 call readi(mcmcard,'N8',n8,0)
1609 call readi(mcmcard,'N9',n9,0)
1610 call readi(mcmcard,'N14',n14,0)
1611 call readi(mcmcard,'N15',n15,0)
1612 call readi(mcmcard,'N16',n16,0)
1613 call readi(mcmcard,'N17',n17,0)
1614 call readi(mcmcard,'N18',n18,0)
1616 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1618 call readi(mcmcard,'NDIFF',ndiff,2)
1619 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1620 call readi(mcmcard,'IS1',is1,1)
1621 call readi(mcmcard,'IS2',is2,8)
1622 call readi(mcmcard,'NRAN0',nran0,4)
1623 call readi(mcmcard,'NRAN1',nran1,2)
1624 call readi(mcmcard,'IRR',irr,1)
1625 call readi(mcmcard,'NSEED',nseed,20)
1626 call readi(mcmcard,'NTOTAL',ntotal,10000)
1627 call reada(mcmcard,'CUT1',cut1,2.0d0)
1628 call reada(mcmcard,'CUT2',cut2,5.0d0)
1629 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1630 call readi(mcmcard,'ICMAX',icmax,3)
1631 call readi(mcmcard,'IRESTART',irestart,0)
1632 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1635 call reada(mcmcard,'DELE',dele,20.0d0)
1636 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1637 call readi(mcmcard,'IREF',iref,0)
1638 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1639 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1640 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1641 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1642 write (iout,*) "NCONF_IN",nconf_in
1645 c----------------------------------------------------------------------------
1646 cfmc subroutine mcmfread
1647 cfmc implicit real*8 (a-h,o-z)
1648 cfmc include 'DIMENSIONS'
1649 cfmc include 'COMMON.MCMF'
1650 cfmc include 'COMMON.IOUNITS'
1651 cfmc include 'COMMON.GEO'
1652 cfmc character*80 ucase
1653 cfmc character*620 mcmcard
1654 cfmc call card_concat(mcmcard)
1656 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1657 cfmc write(iout,*)'MAXRANT=',maxrant
1658 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1659 cfmc write(iout,*)'MAXFAM=',maxfam
1660 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1661 cfmc write(iout,*)'NNET1=',nnet1
1662 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1663 cfmc write(iout,*)'NNET2=',nnet2
1664 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1665 cfmc write(iout,*)'NNET3=',nnet3
1666 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1667 cfmc write(iout,*)'ILASTT=',ilastt
1668 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1669 cfmc write(iout,*)'MAXSTR=',maxstr
1670 cfmc maxstr_f=maxstr/maxfam
1671 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1672 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1673 cfmc write(iout,*)'NMCMF=',nmcmf
1674 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1675 cfmc write(iout,*)'IFOCUS=',ifocus
1676 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1677 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1678 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1679 cfmc write(iout,*)'INTPRT=',intprt
1680 cfmc call readi(mcmcard,'IPRT',iprt,100)
1681 cfmc write(iout,*)'IPRT=',iprt
1682 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1683 cfmc write(iout,*)'IMAXTR=',imaxtr
1684 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1685 cfmc write(iout,*)'MAXEVEN=',maxeven
1686 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1687 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1688 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1689 cfmc write(iout,*)'INIMIN=',inimin
1690 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1691 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1692 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1693 cfmc write(iout,*)'NTHREAD=',nthread
1694 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1695 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1696 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1697 cfmc write(iout,*)'MAXPERT=',maxpert
1698 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1699 cfmc write(iout,*)'IRMSD=',irmsd
1700 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1701 cfmc write(iout,*)'DENEMIN=',denemin
1702 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1703 cfmc write(iout,*)'RCUT1S=',rcut1s
1704 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1705 cfmc write(iout,*)'RCUT1E=',rcut1e
1706 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1707 cfmc write(iout,*)'RCUT2S=',rcut2s
1708 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1709 cfmc write(iout,*)'RCUT2E=',rcut2e
1710 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1711 cfmc write(iout,*)'DPERT1=',d_pert1
1712 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1713 cfmc write(iout,*)'DPERT1A=',d_pert1a
1714 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1715 cfmc write(iout,*)'DPERT2=',d_pert2
1716 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1717 cfmc write(iout,*)'DPERT2A=',d_pert2a
1718 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1719 cfmc write(iout,*)'DPERT2B=',d_pert2b
1720 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1721 cfmc write(iout,*)'DPERT2C=',d_pert2c
1722 cfmc d_pert1=deg2rad*d_pert1
1723 cfmc d_pert1a=deg2rad*d_pert1a
1724 cfmc d_pert2=deg2rad*d_pert2
1725 cfmc d_pert2a=deg2rad*d_pert2a
1726 cfmc d_pert2b=deg2rad*d_pert2b
1727 cfmc d_pert2c=deg2rad*d_pert2c
1728 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1729 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1730 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1731 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1732 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1733 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1734 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1735 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1736 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1737 cfmc write(iout,*)'RCUTINI=',rcutini
1738 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1739 cfmc write(iout,*)'GRAT=',grat
1740 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1741 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1745 c----------------------------------------------------------------------------
1747 implicit real*8 (a-h,o-z)
1748 include 'DIMENSIONS'
1749 include 'COMMON.MCM'
1750 include 'COMMON.MCE'
1751 include 'COMMON.IOUNITS'
1753 character*320 mcmcard
1754 call card_concat(mcmcard)
1755 call readi(mcmcard,'MAXACC',maxacc,100)
1756 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1757 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1758 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1759 call readi(mcmcard,'MAXREPM',maxrepm,200)
1760 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1761 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1762 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1763 call reada(mcmcard,'E_UP',e_up,5.0D0)
1764 call reada(mcmcard,'DELTE',delte,0.1D0)
1765 call readi(mcmcard,'NSWEEP',nsweep,5)
1766 call readi(mcmcard,'NSTEPH',nsteph,0)
1767 call readi(mcmcard,'NSTEPC',nstepc,0)
1768 call reada(mcmcard,'TMIN',tmin,298.0D0)
1769 call reada(mcmcard,'TMAX',tmax,298.0D0)
1770 call readi(mcmcard,'NWINDOW',nwindow,0)
1771 call readi(mcmcard,'PRINT_MC',print_mc,0)
1772 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1773 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1774 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1775 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1776 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1777 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1778 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1779 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1780 if (nwindow.gt.0) then
1781 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1783 winlen(i)=winend(i)-winstart(i)+1
1786 if (tmax.lt.tmin) tmax=tmin
1787 if (tmax.eq.tmin) then
1791 if (nstepc.gt.0 .and. nsteph.gt.0) then
1792 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1793 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1795 C Probabilities of different move types
1796 sumpro_type(0)=0.0D0
1797 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1798 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1799 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1800 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1801 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1802 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1803 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1805 print *,'i',i,' sumprotype',sumpro_type(i)
1806 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1807 print *,'i',i,' sumprotype',sumpro_type(i)
1811 c----------------------------------------------------------------------------
1812 subroutine read_minim
1813 implicit real*8 (a-h,o-z)
1814 include 'DIMENSIONS'
1815 include 'COMMON.MINIM'
1816 include 'COMMON.IOUNITS'
1818 character*320 minimcard
1819 call card_concat(minimcard)
1820 call readi(minimcard,'MAXMIN',maxmin,2000)
1821 call readi(minimcard,'MAXFUN',maxfun,5000)
1822 call readi(minimcard,'MINMIN',minmin,maxmin)
1823 call readi(minimcard,'MINFUN',minfun,maxmin)
1824 call reada(minimcard,'TOLF',tolf,1.0D-2)
1825 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1826 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1827 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1828 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1829 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1830 & 'Options in energy minimization:'
1831 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1832 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1833 & 'MinMin:',MinMin,' MinFun:',MinFun,
1834 & ' TolF:',TolF,' RTolF:',RTolF
1837 c----------------------------------------------------------------------------
1838 subroutine read_angles(kanal,*)
1839 implicit real*8 (a-h,o-z)
1840 include 'DIMENSIONS'
1841 include 'COMMON.GEO'
1842 include 'COMMON.VAR'
1843 include 'COMMON.CHAIN'
1844 include 'COMMON.IOUNITS'
1845 include 'COMMON.CONTROL'
1846 c Read angles from input
1848 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1849 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1850 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1851 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1854 c 9/7/01 avoid 180 deg valence angle
1855 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1857 theta(i)=deg2rad*theta(i)
1858 phi(i)=deg2rad*phi(i)
1859 alph(i)=deg2rad*alph(i)
1860 omeg(i)=deg2rad*omeg(i)
1865 c----------------------------------------------------------------------------
1866 subroutine reada(rekord,lancuch,wartosc,default)
1868 character*(*) rekord,lancuch
1869 double precision wartosc,default
1872 iread=index(rekord,lancuch)
1873 if (iread.eq.0) then
1877 iread=iread+ilen(lancuch)+1
1878 read (rekord(iread:),*,err=10,end=10) wartosc
1883 c----------------------------------------------------------------------------
1884 subroutine readi(rekord,lancuch,wartosc,default)
1886 character*(*) rekord,lancuch
1887 integer wartosc,default
1890 iread=index(rekord,lancuch)
1891 if (iread.eq.0) then
1895 iread=iread+ilen(lancuch)+1
1896 read (rekord(iread:),*,err=10,end=10) wartosc
1901 c----------------------------------------------------------------------------
1902 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1905 integer tablica(dim),default
1906 character*(*) rekord,lancuch
1913 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1914 if (iread.eq.0) return
1915 iread=iread+ilen(lancuch)+1
1916 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1919 c----------------------------------------------------------------------------
1920 subroutine multreada(rekord,lancuch,tablica,dim,default)
1923 double precision tablica(dim),default
1924 character*(*) rekord,lancuch
1931 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1932 if (iread.eq.0) return
1933 iread=iread+ilen(lancuch)+1
1934 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1937 c----------------------------------------------------------------------------
1938 subroutine openunits
1939 implicit real*8 (a-h,o-z)
1940 include 'DIMENSIONS'
1943 character*16 form,nodename
1946 include 'COMMON.SETUP'
1947 include 'COMMON.IOUNITS'
1949 include 'COMMON.CONTROL'
1950 integer lenpre,lenpot,ilen,lentmp
1952 character*3 out1file_text,ucase
1955 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1956 call getenv_loc("PREFIX",prefix)
1958 call getenv_loc("POT",pot)
1959 call getenv_loc("DIRTMP",tmpdir)
1960 call getenv_loc("CURDIR",curdir)
1961 call getenv_loc("OUT1FILE",out1file_text)
1962 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1963 out1file_text=ucase(out1file_text)
1964 if (out1file_text(1:1).eq."Y") then
1967 out1file=fg_rank.gt.0
1972 if (lentmp.gt.0) then
1973 write (*,'(80(1h!))')
1974 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1975 write (*,'(80(1h!))')
1976 write (*,*)"All output files will be on node /tmp directory."
1978 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1979 if (me.eq.king) then
1980 write (*,*) "The master node is ",nodename
1981 else if (fg_rank.eq.0) then
1982 write (*,*) "I am the CG slave node ",nodename
1984 write (*,*) "I am the FG slave node ",nodename
1987 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1988 lenpre = lentmp+lenpre+1
1990 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1991 C Get the names and open the input files
1992 #if defined(WINIFL) || defined(WINPGI)
1993 open(1,file=pref_orig(:ilen(pref_orig))//
1994 & '.inp',status='old',readonly,shared)
1995 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1996 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1997 C Get parameter filenames and open the parameter files.
1998 call getenv_loc('BONDPAR',bondname)
1999 open (ibond,file=bondname,status='old',readonly,shared)
2000 call getenv_loc('THETPAR',thetname)
2001 open (ithep,file=thetname,status='old',readonly,shared)
2002 call getenv_loc('ROTPAR',rotname)
2003 open (irotam,file=rotname,status='old',readonly,shared)
2004 call getenv_loc('TORPAR',torname)
2005 open (itorp,file=torname,status='old',readonly,shared)
2006 call getenv_loc('TORDPAR',tordname)
2007 open (itordp,file=tordname,status='old',readonly,shared)
2008 call getenv_loc('TORKCC',torkccname)
2009 open (itorkcc,file=torkccname,status='old',readonly,shared)
2010 call getenv_loc('THETKCC',thetkccname)
2011 open (ithetkcc,file=thetkccname,status='old',readonly,shared)
2012 call getenv_loc('FOURIER',fouriername)
2013 open (ifourier,file=fouriername,status='old',readonly,shared)
2014 call getenv_loc('ELEPAR',elename)
2015 open (ielep,file=elename,status='old',readonly,shared)
2016 call getenv_loc('SIDEPAR',sidename)
2017 open (isidep,file=sidename,status='old',readonly,shared)
2018 call getenv_loc('LIPTRANPAR',liptranname)
2019 open (iliptranpar,file=liptranname,status='old',readonly,shared)
2020 #elif (defined CRAY) || (defined AIX)
2021 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2023 c print *,"Processor",myrank," opened file 1"
2024 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2025 c print *,"Processor",myrank," opened file 9"
2026 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2027 C Get parameter filenames and open the parameter files.
2028 call getenv_loc('BONDPAR',bondname)
2029 open (ibond,file=bondname,status='old',action='read')
2030 c print *,"Processor",myrank," opened file IBOND"
2031 call getenv_loc('THETPAR',thetname)
2032 open (ithep,file=thetname,status='old',action='read')
2033 c print *,"Processor",myrank," opened file ITHEP"
2034 call getenv_loc('ROTPAR',rotname)
2035 open (irotam,file=rotname,status='old',action='read')
2036 c print *,"Processor",myrank," opened file IROTAM"
2037 call getenv_loc('TORPAR',torname)
2038 open (itorp,file=torname,status='old',action='read')
2039 c print *,"Processor",myrank," opened file ITORP"
2040 call getenv_loc('TORDPAR',tordname)
2041 open (itordp,file=tordname,status='old',action='read')
2042 call getenv_loc('TORKCC',torkccname)
2043 open (itorkcc,file=torkccname,status='old',action='read')
2044 call getenv_loc('THETKCC',thetkccname)
2045 open (ithetkcc,file=thetkccname,status='old',action='read')
2046 c print *,"Processor",myrank," opened file ITORDP"
2047 call getenv_loc('SCCORPAR',sccorname)
2048 open (isccor,file=sccorname,status='old',action='read')
2049 c print *,"Processor",myrank," opened file ISCCOR"
2050 call getenv_loc('FOURIER',fouriername)
2051 open (ifourier,file=fouriername,status='old',action='read')
2052 c print *,"Processor",myrank," opened file IFOURIER"
2053 call getenv_loc('ELEPAR',elename)
2054 open (ielep,file=elename,status='old',action='read')
2055 c print *,"Processor",myrank," opened file IELEP"
2056 call getenv_loc('SIDEPAR',sidename)
2057 open (isidep,file=sidename,status='old',action='read')
2058 call getenv_loc('LIPTRANPAR',liptranname)
2059 open (iliptranpar,file=liptranname,status='old',action='read')
2060 c print *,"Processor",myrank," opened file ISIDEP"
2061 c print *,"Processor",myrank," opened parameter files"
2063 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2064 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2065 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2066 C Get parameter filenames and open the parameter files.
2067 call getenv_loc('BONDPAR',bondname)
2068 open (ibond,file=bondname,status='old')
2069 call getenv_loc('THETPAR',thetname)
2070 open (ithep,file=thetname,status='old')
2071 call getenv_loc('ROTPAR',rotname)
2072 open (irotam,file=rotname,status='old')
2073 call getenv_loc('TORPAR',torname)
2074 open (itorp,file=torname,status='old')
2075 call getenv_loc('TORDPAR',tordname)
2076 open (itordp,file=tordname,status='old')
2077 call getenv_loc('TORKCC',torkccname)
2078 open (itorkcc,file=torkccname,status='old')
2079 call getenv_loc('THETKCC',thetkccname)
2080 open (ithetkcc,file=thetkccname,status='old')
2081 call getenv_loc('SCCORPAR',sccorname)
2082 open (isccor,file=sccorname,status='old')
2083 call getenv_loc('FOURIER',fouriername)
2084 open (ifourier,file=fouriername,status='old')
2085 call getenv_loc('ELEPAR',elename)
2086 open (ielep,file=elename,status='old')
2087 call getenv_loc('SIDEPAR',sidename)
2088 open (isidep,file=sidename,status='old')
2089 call getenv_loc('LIPTRANPAR',liptranname)
2090 open (iliptranpar,file=liptranname,status='old')
2092 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2094 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2095 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2096 C Get parameter filenames and open the parameter files.
2097 call getenv_loc('BONDPAR',bondname)
2098 open (ibond,file=bondname,status='old',readonly)
2099 call getenv_loc('THETPAR',thetname)
2100 open (ithep,file=thetname,status='old',readonly)
2101 call getenv_loc('ROTPAR',rotname)
2102 open (irotam,file=rotname,status='old',readonly)
2103 call getenv_loc('TORPAR',torname)
2104 open (itorp,file=torname,status='old',readonly)
2105 call getenv_loc('TORDPAR',tordname)
2106 open (itordp,file=tordname,status='old',readonly)
2107 call getenv_loc('TORKCC',torkccname)
2108 open (itorkcc,file=torkccname,status='old',readonly)
2109 call getenv_loc('THETKCC',thetkccname)
2110 open (ithetkcc,file=thetkccname,status='old',readonly)
2111 call getenv_loc('SCCORPAR',sccorname)
2112 open (isccor,file=sccorname,status='old',readonly)
2114 call getenv_loc('THETPARPDB',thetname_pdb)
2115 print *,"thetname_pdb ",thetname_pdb
2116 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2117 print *,ithep_pdb," opened"
2119 call getenv_loc('FOURIER',fouriername)
2120 open (ifourier,file=fouriername,status='old',readonly)
2121 call getenv_loc('ELEPAR',elename)
2122 open (ielep,file=elename,status='old',readonly)
2123 call getenv_loc('SIDEPAR',sidename)
2124 open (isidep,file=sidename,status='old',readonly)
2125 call getenv_loc('LIPTRANPAR',liptranname)
2126 open (iliptranpar,file=liptranname,status='old',action='read')
2128 call getenv_loc('ROTPARPDB',rotname_pdb)
2129 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2134 C 8/9/01 In the newest version SCp interaction constants are read from a file
2135 C Use -DOLDSCP to use hard-coded constants instead.
2137 call getenv_loc('SCPPAR',scpname)
2138 #if defined(WINIFL) || defined(WINPGI)
2139 open (iscpp,file=scpname,status='old',readonly,shared)
2140 #elif (defined CRAY) || (defined AIX)
2141 open (iscpp,file=scpname,status='old',action='read')
2143 open (iscpp,file=scpname,status='old')
2145 open (iscpp,file=scpname,status='old',readonly)
2148 call getenv_loc('PATTERN',patname)
2149 #if defined(WINIFL) || defined(WINPGI)
2150 open (icbase,file=patname,status='old',readonly,shared)
2151 #elif (defined CRAY) || (defined AIX)
2152 open (icbase,file=patname,status='old',action='read')
2154 open (icbase,file=patname,status='old')
2156 open (icbase,file=patname,status='old',readonly)
2159 C Open output file only for CG processes
2160 c print *,"Processor",myrank," fg_rank",fg_rank
2161 if (fg_rank.eq.0) then
2163 if (nodes.eq.1) then
2166 npos = dlog10(dfloat(nodes-1))+1
2168 if (npos.lt.3) npos=3
2169 write (liczba,'(i1)') npos
2170 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2172 write (liczba,form) me
2173 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2174 & liczba(:ilen(liczba))
2175 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2177 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2179 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2180 & liczba(:ilen(liczba))//'.mol2'
2181 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2182 & liczba(:ilen(liczba))//'.stat'
2184 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2185 & //liczba(:ilen(liczba))//'.stat')
2186 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2189 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2190 & liczba(:ilen(liczba))//'.const'
2195 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2196 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2197 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2198 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2199 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2201 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2203 rest2name=prefix(:ilen(prefix))//'.rst'
2205 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2208 #if defined(AIX) || defined(PGI)
2209 if (me.eq.king .or. .not. out1file)
2210 & open(iout,file=outname,status='unknown')
2212 if (fg_rank.gt.0) then
2213 write (liczba,'(i3.3)') myrank/nfgtasks
2214 write (ll,'(bz,i3.3)') fg_rank
2215 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2220 open(igeom,file=intname,status='unknown',position='append')
2221 open(ipdb,file=pdbname,status='unknown')
2222 open(imol2,file=mol2name,status='unknown')
2223 open(istat,file=statname,status='unknown',position='append')
2225 c1out open(iout,file=outname,status='unknown')
2228 if (me.eq.king .or. .not.out1file)
2229 & open(iout,file=outname,status='unknown')
2231 if (fg_rank.gt.0) then
2232 write (liczba,'(i3.3)') myrank/nfgtasks
2233 write (ll,'(bz,i3.3)') fg_rank
2234 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2239 open(igeom,file=intname,status='unknown',access='append')
2240 open(ipdb,file=pdbname,status='unknown')
2241 open(imol2,file=mol2name,status='unknown')
2242 open(istat,file=statname,status='unknown',access='append')
2244 c1out open(iout,file=outname,status='unknown')
2247 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2248 csa_seed=prefix(:lenpre)//'.CSA.seed'
2249 csa_history=prefix(:lenpre)//'.CSA.history'
2250 csa_bank=prefix(:lenpre)//'.CSA.bank'
2251 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2252 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2253 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2254 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2255 csa_int=prefix(:lenpre)//'.int'
2256 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2257 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2258 csa_in=prefix(:lenpre)//'.CSA.in'
2259 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2262 write (iout,'(80(1h-))')
2263 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2264 write (iout,'(80(1h-))')
2265 write (iout,*) "Input file : ",
2266 & pref_orig(:ilen(pref_orig))//'.inp'
2267 write (iout,*) "Output file : ",
2268 & outname(:ilen(outname))
2270 write (iout,*) "Sidechain potential file : ",
2271 & sidename(:ilen(sidename))
2273 write (iout,*) "SCp potential file : ",
2274 & scpname(:ilen(scpname))
2276 write (iout,*) "Electrostatic potential file : ",
2277 & elename(:ilen(elename))
2278 write (iout,*) "Cumulant coefficient file : ",
2279 & fouriername(:ilen(fouriername))
2280 write (iout,*) "Torsional parameter file : ",
2281 & torname(:ilen(torname))
2282 write (iout,*) "Double torsional parameter file : ",
2283 & tordname(:ilen(tordname))
2284 write (iout,*) "SCCOR parameter file : ",
2285 & sccorname(:ilen(sccorname))
2286 write (iout,*) "Bond & inertia constant file : ",
2287 & bondname(:ilen(bondname))
2288 write (iout,*) "Bending parameter file : ",
2289 & thetname(:ilen(thetname))
2290 write (iout,*) "Rotamer parameter file : ",
2291 & rotname(:ilen(rotname))
2292 write (iout,*) "Thetpdb parameter file : ",
2293 & thetname_pdb(:ilen(thetname_pdb))
2294 write (iout,*) "Threading database : ",
2295 & patname(:ilen(patname))
2297 &write (iout,*)" DIRTMP : ",
2299 write (iout,'(80(1h-))')
2303 c----------------------------------------------------------------------------
2304 subroutine card_concat(card)
2305 implicit real*8 (a-h,o-z)
2306 include 'DIMENSIONS'
2307 include 'COMMON.IOUNITS'
2309 character*80 karta,ucase
2311 read (inp,'(a)') karta
2314 do while (karta(80:80).eq.'&')
2315 card=card(:ilen(card)+1)//karta(:79)
2316 read (inp,'(a)') karta
2319 card=card(:ilen(card)+1)//karta
2322 c----------------------------------------------------------------------------------
2324 implicit real*8 (a-h,o-z)
2325 include 'DIMENSIONS'
2326 include 'COMMON.CHAIN'
2327 include 'COMMON.IOUNITS'
2329 open(irest2,file=rest2name,status='unknown')
2330 read(irest2,*) totT,EK,potE,totE,t_bath
2333 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2336 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2339 read (irest2,*) iset
2344 c---------------------------------------------------------------------------------
2345 subroutine read_fragments
2346 implicit real*8 (a-h,o-z)
2347 include 'DIMENSIONS'
2351 include 'COMMON.SETUP'
2352 include 'COMMON.CHAIN'
2353 include 'COMMON.IOUNITS'
2355 include 'COMMON.CONTROL'
2356 read(inp,*) nset,nfrag,npair,nfrag_back
2357 if(me.eq.king.or..not.out1file)
2358 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2359 & " nfrag_back",nfrag_back
2361 read(inp,*) mset(iset)
2363 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2365 if(me.eq.king.or..not.out1file)
2366 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2367 & ifrag(2,i,iset), qinfrag(i,iset)
2370 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2372 if(me.eq.king.or..not.out1file)
2373 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2374 & ipair(2,i,iset), qinpair(i,iset)
2377 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2378 & wfrag_back(3,i,iset),
2379 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2380 if(me.eq.king.or..not.out1file)
2381 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2382 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2387 C---------------------------------------------------------------------------
2388 subroutine read_afminp
2389 implicit real*8 (a-h,o-z)
2390 include 'DIMENSIONS'
2394 include 'COMMON.SETUP'
2395 include 'COMMON.CONTROL'
2396 include 'COMMON.CHAIN'
2397 include 'COMMON.IOUNITS'
2398 include 'COMMON.SBRIDGE'
2399 character*320 afmcard
2401 call card_concat(afmcard)
2402 call readi(afmcard,"BEG",afmbeg,0)
2403 call readi(afmcard,"END",afmend,0)
2404 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2405 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2406 print *,'FORCE=' ,forceAFMconst
2407 CCCC NOW PROPERTIES FOR AFM
2410 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2412 distafminit=dsqrt(distafminit)
2413 print *,'initdist',distafminit
2417 c-------------------------------------------------------------------------------
2418 subroutine read_dist_constr
2419 implicit real*8 (a-h,o-z)
2420 include 'DIMENSIONS'
2424 include 'COMMON.SETUP'
2425 include 'COMMON.CONTROL'
2426 include 'COMMON.CHAIN'
2427 include 'COMMON.IOUNITS'
2428 include 'COMMON.SBRIDGE'
2429 integer ifrag_(2,100),ipair_(2,100)
2430 double precision wfrag_(100),wpair_(100)
2431 character*500 controlcard
2433 write (iout,*) "Calling read_dist_constr"
2434 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2436 call card_concat(controlcard)
2437 call readi(controlcard,"NFRAG",nfrag_,0)
2438 call readi(controlcard,"NPAIR",npair_,0)
2439 call readi(controlcard,"NDIST",ndist_,0)
2440 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2441 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2442 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2443 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2444 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2445 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2446 c write (iout,*) "IFRAG"
2448 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2450 c write (iout,*) "IPAIR"
2452 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2456 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2457 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2458 & ifrag_(2,i)=nstart_sup+nsup-1
2459 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2461 if (wfrag_(i).gt.0.0d0) then
2462 do j=ifrag_(1,i),ifrag_(2,i)-1
2463 do k=j+1,ifrag_(2,i)
2464 c write (iout,*) "j",j," k",k
2466 if (constr_dist.eq.1) then
2471 forcon(nhpb)=wfrag_(i)
2472 else if (constr_dist.eq.2) then
2473 if (ddjk.le.dist_cut) then
2478 forcon(nhpb)=wfrag_(i)
2485 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2488 if (.not.out1file .or. me.eq.king)
2489 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2490 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2492 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2493 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2500 if (wpair_(i).gt.0.0d0) then
2508 do j=ifrag_(1,ii),ifrag_(2,ii)
2509 do k=ifrag_(1,jj),ifrag_(2,jj)
2513 forcon(nhpb)=wpair_(i)
2514 dhpb(nhpb)=dist(j,k)
2516 if (.not.out1file .or. me.eq.king)
2517 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2518 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2520 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2521 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2529 if (constr_dist.eq.11) then
2530 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2531 & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
2532 fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2535 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2536 & ibecarb(i),forcon(nhpb+1)
2538 if (forcon(nhpb+1).gt.0.0d0) then
2540 if (ibecarb(i).gt.0) then
2541 ihpb(i)=ihpb(i)+nres
2542 jhpb(i)=jhpb(i)+nres
2544 if (dhpb(nhpb).eq.0.0d0)
2545 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2547 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2548 C if (forcon(nhpb+1).gt.0.0d0) then
2550 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2552 if (.not.out1file .or. me.eq.king)
2553 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2554 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2556 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2557 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2564 c-------------------------------------------------------------------------------
2566 subroutine flush(iu)
2571 subroutine flush(iu)
2576 c------------------------------------------------------------------------------
2577 subroutine copy_to_tmp(source)
2578 include "DIMENSIONS"
2579 include "COMMON.IOUNITS"
2580 character*(*) source
2581 character* 256 tmpfile
2585 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2586 inquire(file=tmpfile,exist=ex)
2588 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2589 & " to temporary directory..."
2590 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2591 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2595 c------------------------------------------------------------------------------
2596 subroutine move_from_tmp(source)
2597 include "DIMENSIONS"
2598 include "COMMON.IOUNITS"
2599 character*(*) source
2602 write (*,*) "Moving ",source(:ilen(source)),
2603 & " from temporary directory to working directory"
2604 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2605 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2608 c------------------------------------------------------------------------------
2609 subroutine random_init(seed)
2611 C Initialize random number generator
2613 implicit real*8 (a-h,o-z)
2614 include 'DIMENSIONS'
2617 logical OKRandom, prng_restart
2619 integer iseed_array(4)
2621 include 'COMMON.IOUNITS'
2622 include 'COMMON.TIME1'
2623 include 'COMMON.THREAD'
2624 include 'COMMON.SBRIDGE'
2625 include 'COMMON.CONTROL'
2626 include 'COMMON.MCM'
2627 include 'COMMON.MAP'
2628 include 'COMMON.HEADER'
2629 include 'COMMON.CSA'
2630 include 'COMMON.CHAIN'
2631 include 'COMMON.MUCA'
2633 include 'COMMON.FFIELD'
2634 include 'COMMON.SETUP'
2635 iseed=-dint(dabs(seed))
2636 if (iseed.eq.0) then
2637 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2638 & 'Random seed undefined. The program will stop.'
2639 write (*,'(/80(1h*)/20x,a/80(1h*))')
2640 & 'Random seed undefined. The program will stop.'
2642 call mpi_finalize(mpi_comm_world,ierr)
2644 stop 'Bad random seed.'
2647 if (fg_rank.eq.0) then
2651 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2652 OKRandom = prng_restart(me,iseed)
2655 tmp=65536.0d0**(4-i)
2656 iseed_array(i) = dint(seed/tmp)
2657 seed=seed-iseed_array(i)*tmp
2660 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2661 & (iseed_array(i),i=1,4)
2662 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2663 & (iseed_array(i),i=1,4)
2664 OKRandom = prng_restart(me,iseed_array)
2667 c r1 = prng_next(me)
2668 r1=ran_number(0.0D0,1.0D0)
2670 & write (iout,*) 'ran_num',r1
2671 if (r1.lt.0.0d0) OKRandom=.false.
2673 if (.not.OKRandom) then
2674 write (iout,*) 'PRNG IS NOT WORKING!!!'
2675 print *,'PRNG IS NOT WORKING!!!'
2678 call mpi_abort(mpi_comm_world,error_msg,ierr)
2681 write (iout,*) 'too many processors for parallel prng'
2682 write (*,*) 'too many processors for parallel prng'
2690 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)